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DETERMINATION OF STABILITY BY LINEAR 
APPROXIMATION OF A PERIODIC SOLUTION 
OF A SYSTEM OF DIFFERENTIAL EQUATIONS 
WITH DISCONTINUOUS RIGHT-HAND SIDES? 
By M. A. AIZERMAN and F. R. GANTMACHER 
(Moscow) 
[Received 28 November 1957] 


SUMMARY 


In the works of Liapounoff (1) the problem of stability of a periodic solution 


i(¢) (With period 7) of a system of differential equations 


dz 
= SilZrsvcepZeqrt) [Lf glSqysvep Sq b+) = Ji(Zqo00-0Z qs b)) (O.1) 
was reduced in the principal (non-critical) cases to the problem of stability of the 
zero-solution a 0 of the system of linear approximation} 
dx “* (cf, 
- Ss { ‘) Xp. (0.2) 
dt hoy OA! z= 0) 
A 


This reduction is applicable when the right-hand sides of equations (0.1) are 
analytic functions, or functions continuous in the neighbourhood of the integral 
curve and differentiable with respect to 2,..., Z,, uniformly in regard to ¢t at the 
points of this curve. 

It is the purpose of this paper to analyse the stability of a periodic solution 
z z(t) t+) 3,(t)] of system (0.1) when the right-hand sides f; are discon- 
tinuous at some points of the integral curve 2z; Z(t), An explanation of the 
meaning of the linear approximation for this ‘discontinuous’ case is offered, and 
theorems, analogous to those of Liapounoff on stability by linear approximation 
in the continuous case, are presented. 


1. Formulation of the problem and principal theorems 

Ix order to obtain conditions to which the right-hand sides of the system 
of equations (0.1) should conform, consider in an (n+ 1)-dimensional space 
Zs.) Z,, £. a curvilinear cylinder C whose axis is the integral curve of the 
undisturbed motion z; = 2,(t). Let an infinite sequence of hypersurfaces 
(discontinuity surfaces)§ 


FAS 5n005 Sut) 0 (1.1) 
+ This article has been published in Russian in Prikl. Mat. Mek. 21 (1957). 
t Here and elsewhere latin indices 7, 7, k have values 1, 2,..., n, and index a values 


§ Henceforth, instead of ‘hypersurface’ and ‘hyperplane’ we will write ‘surface’ and 
‘plane’, 


(Quart. Journ. Mech. and Applied Math., Vol. XI, Pt. 4, 1958) 
S092 .44 ce 
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SUMMARY 
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a =f i(ZyyeeeyZqst) Of g(Zqyeeey Spb +T) = Si(2Zq5---9 Z qs b)] (0.1) 


was reduced in the principal (non-critical) cases to the problem of stability of the 
zero-solution x; = 0 of the system of linear approximationt 


oi - > (2) te (0.2) 


This reduction is applicable when the right-hand sides of equations (0.1) are 
analytic functions, or functions continuous in the neighbourhood of the integral 
curve and differentiable with respect to 2,,.... Z,, uniformly in regard to ¢ at the 
points of this curve. 

It is the purpose of this paper to analyse the stability of a periodic solution 
2, = Z(t) [2,(t+7) = 2,(t)] of system (0.1) when the right-hand sides f; are discon- 
tinuous at some points of the integral curve z; = z,(t). An explanation of the 
meaning of the linear approximation for this ‘discontinuous’ case is offered, and 
theorems, analogous to those of Liapounoff on stability by linear approximation 
in the continuous case, are presented. 


1. Formulation of the problem and principal theorems 


In order to obtain conditions to which the right-hand sides of the system 
of equations (0.1) should conform, consider in an (n+ 1)-dimensional space 
Z,,+++) 2m, ¢ @ curvilinear cylinder C whose axis is the integral curve of the 
undisturbed motion z, = Z,(t). Let an infinite sequence of hypersurfaces 
(discontinuity surfaces)§ 
F,(21).--52,)t) = 0 (1.1) 
+ This article has been published in Russian in Prikl. Mat. Mek. 21 (1957). 


t Here and elsewhere latin indices i, 7, k have values 1, 2,..., n, and index a values 
i am 

§ Henceforth, instead of ‘hypersurface’ and ‘hyperplane’ we will write ‘surface’ and 
‘plane’, 

(Quart. Journ. Mech. and Applied Math., Vol. XI, Pt. 4, 1958) 

5002.44 ce 
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divide the cylinder C into domains H,, and let the integral curve z; = 2,(t) 
intersect the discontinuity surface (1.1) at points M, when t = t,, crossing 
over from the ‘negative’ to the ‘positive’ face of this surface (Fig. 1). 

It is supposed that the right-hand sides of equations (0.1) conform to 
the following conditions: 

(1) The functions f; are periodic in respect to t, with a period r. 

(2) The functions f; are continuous in each domain H,, including the 


,o 














boundaries F,_, = 0 and F, = 0, and on passing through the surface of 
F,, = 0, they have only discontinuities of the first type. 

(3) The functions f; in the domain H, are differentiable with respect to 
Z,,., Z, at the points of the integral curve z; = 2,(t), and uniformly 
differentiable with respect to f¢, i.e. 


SilZ15-++5 ns t) = filE 5 Ft) +> 24 (z;—2;)+0(p), (1.2) 
j 2 


025} 2-20) 


where p= [5 e,—27]' and o(p) >0, asp 0, 
j Pp 


uniformly in regard to ¢t in each interval t,_, < t < #,. 

(4) The conditions of existence and uniqueness of the solution of system 
(0.1) for given initial conditions are fulfilled in every domain H,, and the 
integral curves depend continuously on the initial conditions. As usual, 
it is supposed that the integral curves of the system (0.1) are continuous 
all over C, including the points of intersection with the discontinuity 
surfaces of the right-hand sides f,. 
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As regards the discontinuity surfaces, we make the following assump- 
tions: 

(5) The surfaces F, = 0 are continuous, and at points M, are smooth. 

(6) The derivatives (dF,/dt)- and (dF,/dt)+ along the integral curve 
2; = 2,(t) are non-zero and are of the same sign. 

(7) Translation of the family of surfaces (1.1) along the axis over a 
distance equal to period 7 transfers it into itself. 

Along with system (0.1) consider the system of linear differential 


equations ‘i of 
S on wt 
a <5 >. (2) a (1.3) 
k 


and assume that functions 2; in each interval t,_,+0<t<t,—0 are 
continuous and satisfy equations (1.3) and at t = t, have discontinuities 
determined by the formulae 


ajt—az = €; p32 hy te (1.4.1) 
or aj —az = €; D2 hjay, (1.4.2) 


where €; is the value of the jump of function f; at point M, when crossing 
from the negative to the positive side of surface F, = 0 andt{ 


[e/a ly. 


We will prove in the next paragraph the equivalence of conditions 
(1.4.1) and (1.4.2). These conditions will be called the discontinuity 
conditions. 

The linear differential system determined by the linear differential 
equations (1.3) and by the linear discontinuity conditions (1.4) will be 
called the linear approximation of the initial non-linear system (0.1). 

The significance of the linear approximation thus introduced will be 
seen from the following theorem. 


THEOREM 1. If the zero-solution of a system of linear approximation is 
asymptotically stable, then a periodic solution z, = Z,(t) of the initial non- 
linear system (0.1) is also asymptotically stable. 


The characteristic equation of linear approximation is defined in the 
same way as in the continuous case. Denote the fundamental matrix of 





+ Here os “(3% 2 f,+ Fe 


ze s=4it)” 
and indices — and + define values taken at t = t,—0 pe t = t,+0, i.e. before and after 
the point of discontinuity. 
t Quantities ¢; and hj depend also on index a, but for the sake of clearness it is omitted. 
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the linear approximation by X(t). Its columns are linearly-independent 
solutions of the systems (1.3)-+-(1.4). Then the following identity holds: 


X(t-+7) = X()U, (1.6) 
where U is a certain non-singular constant matrix. The characteristic 


equation det(U—pH) = 0 of the matrix U (here £Z is a unit matrix) is 
called the characteristic equation of linear approximation. 


THEOREM 2. If the absolute value of at least one of the roots of the 
characteristic equation of linear approximation is larger than unity, then the 
periodic solution z; = 2,(t) of system (0.1) is unstable. 


2. Preliminary remarks 


1°. The peculiarity of Theorems 1 and 2 lies in the fact that the analysis 
of the stability of continuous trajectories of system (0.1) is replaced by 
the analysis of the stability of discontinuous trajectories of the linear 
approximation (1.3)+-(1.4). It is assumed that for these discontinuous 
trajectories the definitions of stability and asymptotic stability are the 
same as those introduced by Liapounoff for continuous trajectories. 


2°. All the equations cited above were written in scalar form in order 
to be able to formulate the generalization of Liapounoff’s theorem, in the 
same terms and notation as in Liapounoff (1). In order to prove these 
theorems it is more convenient to state the basic relations in matrix form. 
For this purpose consider columns z, f, x, &, h-, and h+, composed of 
elements z;, f;, 2;, €;, hj , and hj respectively, and also square matricest 


d mal B- = |é,hzl| = gh", BY = [hg || = gh’. 


Now we can write the system of initial equations (0.1) and the linear 
approximation (1.3)+-(1.4) in matrix form: 


d 


—, = Set) (2.1) 


dx of 
nee ee Ee 3 2.2 
dt (z ) a , 


at+—xz- = B-x- (2.3.1) 


or at+—xz- = Bret. (2.3.2) 


3°. Let us use the matrix form of notation in order to prove the 
equivalence of the discontinuity conditions (2.3.1) and (2.3.2), or (1.4.1) 


+ Here and in what follows the accent ’ denotes the transposed matrix. 
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and (1.4.2), which are the same. Equalities (2.3.1) and (2.3.2) may be 
written in the form: 


a+ = (E+B-)cr-, a+ = (E—Bt)-12-. (2.3.3) 
We will show that 2+ B- = (E—B*)-1. Indeed, this equality can be 
reduced to equality B-— B+— B B+ = 0, the correctness of which follows 
am B-Bt = th~'th’ = (hh = cBY (c = h~“€) 
and B-—Bt = | aan ai ar ar = cBe. 
Here éF,/éz, i.e. the column of elements @F,/0z,, is the gradient of function 


F,. From the equivalence of conditions (1.4.1) and (1.4.2) it follows that 
whenever ¢ = t,, we have 


p2 hz x, = p32 hj af. (2.4) 
Now consider the matrix 


S = E+ B- = (E—B+)-. 
We will write the conditions of discontinuity in the form 


at = Sz. (2.5) 


Let us now examine the structure of matrix S and the transformation (2.5) that 
it generates. Let K denote the intersection of the plane tangential to the surface 
of discontinuity F, = 0 (at the point M,) with the plane t = t,. The vectors in 
the plane ¢ = ¢,, parallel to K satisfy the equation h~’x = 0. 

_ (GP, /dt)ing 
¥* GP ding 
In accordance with condition (6) of section 1, the quantity y > 0. 

Let us examine the three cases separately. 


Case 1. y #1. In this case ¢c = h~’§ # 0 and £ does not belong to K. For any 
vector x we have a unique expansion 2 = xg+2, (x¢||&, 2, € K). 
Since SE = £+£(h~’€) = y& and x¢\|£ we have 
Sig = ye. 
On the other hand, Sa, = 2% +&(h~’xy) = Ly. 


Let = 1+h~€ = l+e. 


Hence, the transformation (2.5) does not change vectors parallel to K, but vectors 
parallel to € are multiplied by y. Between the components of vectors z~ and 2+ 
there exists a relation + + + re 

at = yzz, af = 2. 


It follows that matrix S has simple elementary divisors and the characteristic root 
y is simple and | is of (n—1)-multiplicity. 

Case 2. y = 1, #0,h #0. Here €e K. All characteristic roots of the matrix 
S are equal to 1, but S # HZ. Matrix S has one elementary divisor of the second 
degree, all the others being of the first. 

Case 3. y = 1 and = 0 orh = 0, In this case S = E and z+ = z-. Therefore 
all integral curves of the linear approximation are continuous at * = t, only when 
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the right-hand sides of the differential equations (0.1) are continuous at point M, of 
the integral curve z; = 2,(t), or if at this point the discontinuity surface is tangential 
to plane t = ?,. 


4°. Let us introduce the discontinuous functions 
oF ,/éz 
4,(t) = |= ; 
ul , coca 
It follows from (2.4) that the function D3 h,(t)x,(t) is continuous. There- 


fore the system of linear approximation may be written 
Of; > 
Dx, => (2) +36 in(t) (tt) | (2.6) 
TUNE s=) = 
where 8(¢) is the Dirac-function and Dx, the derivative of function x; with 
regard to 5-terms.t 
On the other hand, it can be shown by simple calculations, that the 
expression in square brackets is a generalized derivative with respect to 
z, of the right-hand side f; (i.e. a derivative with regard to 5-terms) 
(Defdewo = (2) + ¥ &hy(t)3(t—t,). (2.7) 
“k/ 2=2) a 
Thus, the system of linear approximation (1.3)+-(1.4) can be written 
Du, = y 3 (Dy fiens%r (2.8) 


which differs from the linear approximation (0.2) with continuous right- 
hand sides only} in having the ordinary derivatives éf;,/0z, replaced by 
generalized derivatives D,f;, and the product of a continuous function 
u(t) by 5(t—t,) is equal to u(t,)d(t—t,). It should be remembered that 
€, and h,(t) depend on index «, which, as has already been pointed out, is 
omitted. 
5°. We will now consider the deviation 

r= z—2(t) (2.9) 

and rewrite equation (2.1) in terms of deviations 


- = p(x, t) [ p(x, t) = f(2(t)+-a, t)—f (X(t), t)]. (2.10) 


+ We proceed here from the fact that the derivative of any function g(t) that is dis- 
continuous at ¢ = t, with a jump 7, may be expressed with regard to 5-terms by a sum 


dg 
—+¥F t— . 
at Na d(t—t,) 


{ As applied to a particular case of a system, differing from a linear system with constant 
coefficients by having on the right-hand side one non-linear function of one argument of 
the type sign z, in (2), the linear approximation was composed of 5-functions in a form 
which is a particular case of equation (2.6). The author based his work, in essence, on the 
theorem of Liapounoff, although the application of it in the case of discontinuous systems 
was not proved. 
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The equations of the discontinuity surfaces will now be 
Pia,t)=0 [P,(z,t) = F(Xt)+2,1)]. (2.11) 
In (n+1)-dimensional space (x,t) the discontinuity surface (2.11) is 
continuous, and so is the surface F, = 0 in the space (z,¢), but in contrast 
with the surface F, = 0 in the space (z,t), the surface (2.11) in the space 
(x,t) is not smooth and has an edge at the line of intersection with plane 
t = t,. 
Because of the second term in formula (2.10) the function p(z,t) has 


a" 


| 
| Mas ji! 
| tact! Pees 


/ | 
Ma ff! 


S= 


t Ene T 
l 
| 
| 
| 








/ 


M,sif7/ 
als,’ 





| 
/ | Pa 
| 


te-1 ; 





- / Pp 
» Trajectory of | | a-1 
disturbed motion Ls Trajectory of. ‘ 
7 linear approximation 








Fie. 2 


a discontinuity not only on surfaces P, = 0, but also on planes t = t,. 
The solutions of systems (2.10) are also continuous, as are the solutions 
of the initial system (0.1). The space (z,t) is transformed into the space 
(x,t) (Fig. 2) and each plane ¢ = const. is translated upon itself. The point 
of intersection of this plane with the curve z; = 2,(t) transforms into a 
point of the t-axis. Now, the t-axis becomes the axis of cylinder C in 
space (x,t). The surfaces P, = 0 and planes t = t, divide the cylinder C 
into central and angular domains. The latter are located between the 
surfaces P, = 0 and the corresponding planes t = t,. We will name the 
angular domains upper or lower, depending upon whether the boundary 
surface P, = 0 of the domains lies over or under the corresponding plane 
t = t,, respectively. At each point of the t-axis, differing from points 
M,, the function p(z,t) is zero. At the point M, of the t-axis (¢ = t,) the 
function p(x,t) may have three values: zero when approaching from 
the central region, € and —é when approaching, respectively, from the 
lower or upper angular domains. 

6°. Now we can give a simple geometrical interpretation of the dis- 
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continuity conditions. Let us suppose, for instance, that an integral curve 
of the non-linear system (2.10) intersects the lower angular domain, 
adjoining the plane t = t,. We will choose the trajectory close to the 
t-axis. Thus, we may consider that its direction is nearly ‘vertical’. Let 
us continue the trajectory through an angular domain, maintaining its 
‘verticality’. Denote the point of its intersection with the plane ¢ = t, 
by x~ and the point of intersection of the same plane with a true (broken) 
integral curve by 2+, then, to within infinitesimals of higher order, we 
will obtain from the right-angled triangle the discontinuity conditions (2.5). 

Indeed, let us replace the equation of the discontinuity surface P, = 0 
by the equation of the tangential plane 


(Fe) 2+ (Ge) t-te) = 0. 


0x } 4, at }u, 


From (2.11), it follows that 


OP)’ _ (0s oP, _ (dF) 
Ott, \ a Iau. et Jy, \ at ay 


By substituting these derivatives into the equation obtained for the 
tangent plane, and solving it for t—t,, we get 





t,—t = h-'x. 
Ignoring infinitesimals of second order, the inclination of the trajectory 
in the angular domain is £. Thus 
at—a~ = §(t,—t,), 
where ¢, is the value of ¢ at the entrance of trajectory into the angular 
domain, i.e. at x = x-. Hence 
at—2z- = th-'z- = B-z-, 


thus proving the point. 


3. The Liapounoff transformation 
Apply a Liapounoff transformation 

a = L(t)y, (3.1) 

where L(t) = X(t)e-“*, A= “In U, (3.2) 


to the linear approximation. Here L is the transformation matrix, X the 
fundamental matrix of the linear approximation (2.2)+(2.3), U the 
constant matrix in identity (1.6), y the column, composed of new variables 
Vay-orr Fat 


+ Transformation (3.1) is called a Liapounoff transformation if all the elements of 
matrices L, L~, and dL/dt are bounded, when t > fp. 
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Let Y denote the fundamental matrix of system (2.2)+-(2.3), transformed 
to variables y;. It then follows from (3.1) that X(t) = L(t)Y(t), or, taking 
into account (3.2), we get Y(t) = e4*. 

It follows directly from this that the equation of linear approxima- 
tion in variables y; becomes 


— = Ay (3.3) 


and the integral curves of the system of the linear approximation in the 
space (y,¢) are continuous. Previously, the linear approximation in space 

















t \ 
at 
\ 
‘\ 
| \ t=ba+t 
tari 1 Qa 
\ 
\ t=te 
ta 1 
Qa 
gg 1 
\ 
\ t=ta-1 
ta-t ‘ 
Y a 
Trajectory of a-! 
disturbed motion i Trajectory ” Ss 
// inear approximation 
| _ 
Fie. 3 


(x,t) was described by the linear equations (2.2), with variable piecewise 
continuous coefficients and by the jump conditions. Now in the space 
(y,t) the jump conditions are absent, because of the continuity of the 
integral curves. 

Thus, the linear approximation investigated in this paper may be 
transformed to a system with constant matrix of coefficients; but unlike 
the case dealt with by Liapounoff, the periodic matrix L(t) is discontinuous 
at t = t,, since at these values of ¢, the columns of matrix X(t) are dis- 
continuous. Indeed, by (2.5) X+ = SX-, when t = t,, and hence 


L+ = SL-. (3.4) 
Applying the Liapounoff transformation to the initial non-linear system 
(2.10), we obtain the system 


dy _ a» tia _4L 


The discontinuity surfaces P,(z,t) = 0 are transformed into surfaces 
Q.(y,t) = 0, where Q.(y, t) = Pi(L(ty, t) (Fig. 3). 
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Unlike the surfaces P,(x,t) = 0, which are continuous and on the lines 
of intersection with the planes t = t, having edges only, the surface 
Q.,(y,t) = 0 has discontinuities along these lines of intersection, i.e. it is 
formed out of two continuous parts (at ¢ < t, and ¢t > t,) intersecting the 
plane ¢t = t, along different curves, which pass through one point t = t, 
of the t-axis (corresponding to point M,). In the coordinates (x,t) the 
integral curves of the system (2.10) are continuous. The integral curves 
of (3.5), corresponding to these curves in space (y,t) are discontinuous at 
t = t,, since at these values of ¢ the matrix L(t) is discontinuous. Att = t, 
we have 


= [tyt = Ly; (3.6) 
hence, on account of (3.4) 


y* = (L-)1"8S1L-y-. (3.7) 


Thus, in the space (y,t) the integral curve of undisturbed motion, as in 
space (x,t), coincides with the t-axis. The curve of disturbed motion in 
the space (y,?¢), as in space (x,t), has edges (on the surfaces Q, = 0), but 
in the space (y,¢), unlike the space (x,t), it also has discontinuities on 
planes t = t,, defined by formula (3.7). The trajectories of the linear 
approximation in the space (y, t) are continuous and are defined by a system 
of linear differential equations (3.3) with coefficients constant everywhere. 
The space (y,t), as also the space (x,t), is divided by the surfaces Q, = 0 
and the planes ¢ = ¢t, into central and angular domains. 


4. Proof of the theorems 


Let the zero-solution of the linear approximation (1.3)+-(1.4) be 
asymptotically stable. Then the solution y = 0 of system (3.3) is also 
asymptotically stable, and all characteristic roots of matrix A have 
negative real parts. In this case, as was shown by Liapounoff, there exists 
a positive-definite quadratic form V(y) whose total derivative (dV /dt)®, 
calculated by means of (3.3), is negative-definite: 

, dv\° 
Vy) > 9, (=) <0 (y¥ #9). 
dt 

Let us consider the variation of values of this function V(y) along 
discontinuous integral curves of the non-linear system (3.5). 

An integral curve passes through central and angular domains and has 
discontinuities on the plane t = ¢t,. We will now show that the value of 
function V along trajectories in the central domain decreases and we will 
estimate the rate of decrease. In angular domains the value of V may 
increase. We will estimate the rate of increase and show that the increase 
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is compensated by jumps of function V on planes t = t,, adjoining angular 
domains. 

In the central domain the linear approximation of non-linear system 
(3.5) is the system (3.3). Thus, in the central domain 


q = Ay+(*). (4.1) 
Here and subsequently (*) and (**) stand for the sum of terms, the order 
of smallness of which, in regard to y is greater than one or two respectively. 


Then 
1/av\ 1/dv\? 1 
ra) = 7a) +H") si 


where derivatives (dV /dt) and (dV /dt)® are calculated in accordance with 
systems (3.5) and (3.3) respectively. 

The first term of the right-hand side of (4.2) is a continuous, homo- 
geneous function of zero degree and it has a negative maximum —p?. 
The second term tends to zero when y > 0 and therefore, for any p? < pu? 
and y sufficiently small, in consequence of (4.2), 


1 /dv\" 
wa) < —p}, (4.3) 
and thence V < V*e-¥it”, (4.4) 


Here V* and V denote the values of V at moments ¢* and t (t* < ¢t) when 
the point of the integral curve is in the given central domain. Formula 
(4.4) defines the rate of decrease of function V along the integral curve in 
the central domain. 

In the angular domain as y - 0 (¢ then also tends to the respective t,), 


the function g has finite limits y~ = (L-)“€ and n+ = —(L*)~€ in the 
lower and upper angular domains, Se, Thus, in an angular 
domain 1 dV 


eV 
Vid Vi nD by = Pi ne 
is bounded, when y is small enough, i.e. there is a K » ch that 


1 dV 
2K; 

yi dt|~ 
hence \Vi—(V’)t| < Ki\t—t’|, 
where V and V’ are the values of V at t and ¢’ in one and the same angular 
domain. Since |t—t’| is a quantity at least as small as the first order of 
magnitude in respect to y, the ratio V/V’ is bounded when y is small 
enough, i.e. there exists a number N > 1, such that 


V <NV’', (4.5) 
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i.e. the coefficient of increase of V in an angular domain is finite. We will 
now show that the variation of V in an angular domain is compensated 
by a corresponding jump V(y+)—V(y-). Assume, for instance, that the 
integral curve passes through the lower angular domain in spaces (z, t) 
and (y, t) from points 2, t, and y,, t, on the discontinuity surface to points 
(z,t,) and (y~,t,) respectively on plane t = t, (Fig. 4). 
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t 4 
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Fic. 4 


Solve the equation of surface P,(x,t) = 0 with respect to t:+ 
t,—t = h-’x+(*). (4.6) 


Since, in the lower angular domain lim p = é, as x > 0. 
te 
%—2, = [p dt = &(t,—t,)+(*) 
th 


= th~'x,+(*) = B-z,+(*), 
1.€. z= Sx,+(*). (4.7) 
On the other hand, 
= L(ty)y, on L-y,+[L(t,)—L(t,—9) Jy, = L-y,+(*). (4.8) 
From (3.4), (3.6), (4.7), and (4.8) we have 


Ltyt = 2 = Sx,+(*) = SL-y,+(*) = Ly, +(*). (4.9) 
By multiplying from the left by (Z*+)-, we get 
yy =nt(*), Viyt) = Viy,)+(**) (4.10) 


+ Since z = Ly, y = L~'z, and the elements of the matrices L(t) and L~(t) are bounded, 
quantities infinitesimal in respect to y will be also infinitesimal in respect to x, and of the 
same order of magnitude, and vice versa. Hence we may use here the same designations 
(*) and (**). 
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and thus, when y is small enough, 


eS 
Viy) Tw 


where 7» is a small positive quantity of the first order. 

Let 7 = min,(t,,,—t,); let v be an arbitrary positive number, smaller 
than pw, and let »y < p, <p. Assume that 7 < }(y?—v*)T. Let « > 0 be 
so small that inside the elliptic cylinder V = « the inequalities (4.4), (4.5), 
and (4.11) hold good, and the interval of time At spent inside any angular 
domain contained in the cylinder V = « is smaller than }(u?—v?)T'/p?. 

Such a choice of ¢ can be achieved, since it follows, from the periodicity 
and property (7) (section 1), that the number of central and angular 
domains under consideration is finite. Consider a sequence of moments 
of time t® = }(t,+t,,,). The corresponding planes t = ff inside cylinder 
V = e do not cross angular domains. 

Let 5 = «/N <e. Take an initial point of an integral curve at t = tf 
inside the cylinder V = 8. This point belongs to the central domain. 
When ¢ varies over the interval tf < t < ¢§, function V at first decreases 
along an integral curve; then it may increase (the increase is compensated 
by a jump), and again decrease. Since the coefficient of increase does not 
exceed NV, the integral curve on segment tf < t < ¢ lies inside the cylinder 
V = «. 


Now, assuming that V, = V._., we have 





ct < (**) < e”, (4.11) 


Ve < Viewkr-son < Vet. (4.12) 


Since, V, < V, < 3, it follows that the point of an integral curve at 
t = tf will again be inside the cylinder V = 5, and thence we can repeat 
the foregoing argument for the next interval t{ < t < ¢} of the integral 
curve, and so on. It follows, that any integral curve of the non-linear 
system, beginning at t = tf inside cylinder V = 3, will always be inside 
cylinder V = «, and Vz < Vie-o-r, 


Thus, lim V, = 0 when «ao, Since the inequality V < NV, holds good 
in every interval  <t < f,,, lim V(t) = 0 when to. Thus, Theorem 
1 is established. 

Starting to prove Theorem 2, we assume that the absolute value of at 
least one of the roots of the characteristic equation of the linear approxima- 
tion is greater than unity. Then one of the roots of the characteristic 
equation of matrix A has a positive real part. In this case there exists 
a quadratic form V(y) which is positive for a certain value y, and, by (3.3), 
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has a positive-definite derivative 


. , dv\° 
h=Vy)>% (7) >0 #9). 
' . 1/dV\° 
Let «x? = min va) (V > h). 
Then, in any central domain along an interval of an integral curve with 
V > \ we have 
1 (dV\" 
ain. Seal Ali 
Aa) =— i, (4 13) 


when y is small enough. 

Here 0 < x, < x, and (dV/dt)™ is a derivative, calculated by equation 
(3.5). Again, let 7’ = min,(t,,,—t,) and 0<o@ < x?. 

Let « > 0 be small enough, so that inside the cylinder V = « the 
inequalities (4.11) and (4.13) hold good, and let the time At spent in any 
angular domain contained in cylinder V = ¢« be smailer than 
sg (xi—o*)T. 


» 
aK) 


Assume » < 4(«x?—o*)7’. Then, instead of inequality (4.12) we have 
V, > V, et+*KT-4)-1 > V, e?T {V, = V(y)}- 

This inequality will hold good, if at all values of the interval tf < t < #, 
the integral curve lies inside the cylinder V = «. Then V, > V, > 0 and 
we may repeat our argument for the next interval of the integral curve, 
ete. At a certain moment of time the integral curve will leave the cylinder 
V = e, since the values of V, increase quicker than the terms of a geo- 
metrical progression with a common factor g = e*? > 1. As V(y,) can 
be made as small as we like, Theorem 2 is established. 
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SUMMARY 


It is shown that the expansion assumed by Goldstein to describe the flow near 
separation in a laminar boundary layer is incomplete and that further terms which 
include powers of logarithms must be added. These terms are individually singular 
at separation. Although it cannot be inferred that the velocity profile must also be 
singular at separation, it is suggested that if the boundary layer is to continue 
downstream of separation the main stream must adjust itself so that these terms 
cannot appear. The solution may then be continued through separation by means 
of a power series into a region of reversed flow. However, it is shown that in 
addition to the power series an infinity of new terms may appear in the solution 
downstream of separation which is therefore no longer specified uniquely by the 
mainstream velocity and the velocity profile at the beginning of the boundary layer. 


1. Résumé 


Ten years ago Goldstein (1) published an important paper on the flow 
in a two-dimensional boundary layer near separation, i.e. near the point 
where the skin friction vanishes. He considered the non-dimensional form 
of the boundary layer equations 


Cu, Ou _ ytl¥ Ou any __ Op 


ne ONE ee ES SP r 38? u=>=-z-> v=; 
ex Gy dx ey oy ex 


u 





(1) 


in which x measures distance upstream along the wall from the point of 
separation O and y distance normal to the wall as shown in Fig. 1. 
Further, u, v are the velocity components in the directions of x decreasing 
and y increasing respectively and y& is the stream function. The boundary 
conditions are that u = v = 0 at y = 0 and that wu tends to prescribed 
functions U(x), k(y) as yoo and as x tends to some positive value 
respectively. The non-dimensional form chosen was such that 
U(0) = U'(0) =1 
and therefore near x = 0 we may write 
dU : 
U—=1 Pw, 2 
z= 1+ EF (2) 
where P. are prescribed constants. 


(Quart. Journ. Mech. and Applied Math., Vol. XI, Pt. 4, 1958] 











400 K. STEWARTSON 


Goldstein assumed that at x = 0, where the skin friction vanishes and 
the boundary layer separates, 


u= yt Lay, (3) 


where the a, are unknown to begin with. The coefficient of y is zero in (3) 
because (0u/éy),.9 = 0 at separation and the coefficient of y? is 4 as a 


Ula) 


Outer limit 
rod boundary 


layer 3 a 





ooo” 








——Reversed flow 
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Fic. 1. Velocity profiles in the boundary layer near separation 


result of the choice of units. His object was to determine the behaviour 
of % in the neighbourhood of x = 0 by solving (1) with boundary condi- 
tions u = v = 0 at y = 0, and (3). Although thea, are unknown initially, 
if his procedure were correct, they would all be determined in terms of 
the P’s and one other parameter. His method was to write 


E=2t,  g=yl2z), $= PY PF), (4) 


substitute into (1) and equate éaiiape: of £. The boundary conditions on 
f, are that 





f.(0) =f%(0) =0 and lim fe t= a... (5) 
aid 
As a result of the solution he found that 
ou = 
ated =n 9 4(r+1) 6 
Ge. 2” aa 


where a, = }f7(0). 
On substituting (4) into (1) and comparing coefficients of ¢ it was found 


om fo(n) = $7’, Pca | 7*, (7) 
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and f, satisfied 


fr 40°F +h(r+4)nf,— (7 +3)nf, = Gn) (r > 2). (8) 


The right-hand side G, of (8) is a known function of f, when s < r and, 
if }r is an integer, of the pressure coefficient P,,; it has the form 


—2(r+2)ay a, 7?+GE(m), (9) 


where G, is explicitly independent of «,_, and has a double zero at the 
origin. 

Since one of the complementary functions of (8) is y? it is clear that a, 
cannot be determined so long as the expansion of ¥ is not continued 
beyond f,. Goldstein (1), however, suggested a means by which they 
could all be expressed in terms of the P’s and a,, and with Jones (2) 
evaluated a,, a3, and a, as follows. The other two complementary func- 
tions of (8) are g,, h, of which g, has a simple zero at » = 0, h,(0) 4 90. 
Both g,, h, may be exponentially large at infinity but a suitable combinaton 


of them can always be found which is algebraic at infinity. The appro- 
priate solution of (8) is 


f, = 40, ~,4(n—9,) +5), 


where j, satisfies (8) with G, replaced by G, and hence is independent of 
a,_,;. Now if we choose j, to have a double zero at 7 = 0 it will in general 
have an exponentially large component at infinity which can be cancelled 
by a suitable multiple of g, provided only that g, is also exponentially 
large, i.e. r ~ 4n+2, where nis an integer. Hence if the solution is known 
as far as f,_, in terms of a, and the P’s, (8) will determine «,_, except if 


r = 4n+-2. In the latter case however g, is a polynomial and instead G, 
must satisfy the condition 


3) 


| G.e-#n'(2ng,—n%g,) dn = 0. (10) 


0 


The first difficulty occurs when r = 6 and Jones (2) showed that in this 
case the value of the integral in (10) is —4a$+4a§ so that it is distinctly 
possible for (10) to hold when r = 6. Goldstein pointed out that in this 
case the way is clear for a complete determination of , for a, would then 
quite possibly be determined from the integral condition at r = 10 and 
so on. Thus all the a, could be found in terms of «, and the P’s, and 
presumably «, could then be found from the condition that u— U as 
yoo. Since f, would be completely determinate all the a, would be 


5002.44 pd 
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known and all separating boundary layers with the same main stream 
would have the same separation profile. 


2. Discussion 


Among the many aspects of this important work there are two of special 
interest here. First, does the condition (10) really hold when r = 6? 
Jones’s work indicates that it may possibly not hold and further it is 
unusual to find a mathematical theory of this kind hanging by a thread. 
Moreover in his corresponding study of laminar wakes, Goldstein (3) 
encountered a similar integral to (10) which he showed to be non- 
vanishing. 

Second, if «, determines the asymptotic behaviour of % near separation 
why is it that when a, = 0 an infinity of arbitrary constants appear in 
the solution (1, section 6), a,,,, then being arbitrary? 

These questions may be answered by making use of a theory of asymp- 
totic expansions given in a previous paper (4). There it is shown that 
certain asymptotic expansions, developed in connexion with problems in 
boundary layer theory, are incomplete and it is also shown how the 
additional terms necessary to complete them may be calculated. Asso- 
ciated with the new terms, however, are arbitrary constants, which cannot 
be determined by the method and depend in some way on the initial 
conditions. 

Before doing this however it is necessary to change the point of view 
slightly. Goldstein’s plan was to assume that w could be expanded as a 
power series in y at x = 0, the point of separation, and deduce y as a 
power series in € = x!, whose coefficients were functions of » = y/2tzt, 
which would be valid in zx > 0, y > 0. Unfortunately, the modification 
to u, which is necessary, contains terms which tend to infinity as « > 0+ 
for fixed y > 0. The initial assumption that wu can be expanded as a 
power series in y when x = 0 is therefore invalid. 

The expansion for y obtained is however not necessarily invalid in the 
entire neighbourhood of the origin. For Goldstein could equally well have 
assumed, instead of (3), the expansion of (@u/éy),.. as the series of 
ascending powers of x! which is given in (6), deducing from it and the 
other boundary conditions the behaviour of % in the neighbourhood of 
the origin. Looked at from this point of view the aim would then be to 
determine if possible the behaviour of u when xz = 0 and y #0. The 
modification to his solution adds more terms to his expansion for (@u/@y),,_9 
without impairing the validity of the basic notion. Since some of the new 
terms tend to infinity as x > 0+ for fixed y > 0, however, it is no longer 
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possible to extend the range of validity of the expansion for y as far as 
this line. It is only valid in fact when 

x>0, y < Az, 
where A is any finite positive number. 

The problem may be posed in the following way. The simplest expan- 
sion of % about the origin is a regular double power series in x and y. 
According to it (@u/éy),.9 may be expanded as a series of integral powers 
of x near x = 0. However, in the only numerical solution of the boundary 
layer equations in the vicinity of separation which is available, (@u/éy),,_9 
is not of this form. How can the expansion of y be generalized to account 
for the computed behaviour of (@u/éy),.9? The correct way is to assume 
that y% can be expanded as a series of powers of € and log whose coeffi- 
cients are functions of » which can be determined, apart from certain 
numerical factors, by substitution in the differential equation (1). 


3. The modification to the theory 
We shall assume as a new starting-point therefore that 

ri: _ oF > phr+y 

fF). ; 2 is °) 
and the appropriate boundary conditions at 7 = 0 are 

fA0) =f(0) = 9, — f(0) = $a, 

In order that the motion at x = 0+-, y # 0 should not be violent we must 
require in addition that f, be not exponentially large at infinity. The a, 
are not now all independent and it follows that if (6) is the correct expan- 
sion for (@u/éy),-9 then the new conditions are equivalent to (5). The 
solution obtained with the new standpoint is thus equivalent to the other 
and leads to the same difficulties. We now ask whether (6) is the most 
general form for (@u/é@y),-9; can more terms be added to (6) without im- 
pairing the validity of the corresponding expansion of 4? Using the 
results in (4) these questions may be answered and further the teory 
rendered self-consistent. These imply that an infinite number of infinite 
series of ascending powers of € whose coefficients are the products of a 
power of log é and a function of 7 can be added to (4) and correspond- 
ing series to (6). Each series is associated with one of the integral con- 
ditions (10) which occurred in Goldstein’s original theory in a way which 
is exemplified by the case r= 6. Knowing that the new terms may 
have logarithmic factors it is found by trial that the most general form 
of % possible is 


$= 216 S Ef(n)-+24EFlog é Fal) +€F(n)]+ 01 log é), (11) 








404 K. STEWARTSON 





whence 


u = 2¢2 > & f'.(n) +2€7 log &[ F's(n) +€F ¢(n)]-+ 0(€* log £) 
and 


(=) = 2 ¥ E°f(0) 4286 log Ef F5(0) +-£F(0)] + O(log é). 
y=0 : 


The reason is as follows. First let the term of lowest order which can be 
added to (4) be &*+*F(n) where s > 0. Then F,(n) = a, y?, where a, is 
a constant. This new term initiates a series of terms just as 24a, £4? does 
in the original series. The equation for the next term in the new series, 
which is proportional to £+® or £4+* according as s < 1 ors > 1, is of the 
form (8) but the right-hand side is proportional to 72. Hence G = 0 and 
the term is exponentially large at 7 = oo except if a, = 0 or ifs = 4r+1, 
r integer, when it is a polynomial, and £*+*F,(n) is already in (4). Thus 
no extra terms proportional to powers of £ can be added to (4). Next let 
the term of the lowest order which can be added to (4) be &+(log €)"F,() 
where t > 0. As before F(n) = «7? where «, is a constant and ¢t = 4r+-1, 
r integer. The constant a, is not now determined by the equation for the 
term proportional to £**(log é)* but by the equation for the term pro- 
portional to €**(log €)"-!, which leads to an integral condition like (10). 
However, the G is zero, except if t, = 1, r == 1, 2,..., when it depends on 
f.{n), 8 < 4r+-1. Thus a, need not be zero if t = 4r+-1, t, = 1, accounting 
for the terms in (11). 

Substituting (11) into (1) we find that f,, r < 6, satisfies the same 
equation as before. The equation for F;, is 


F; —49°F5+37°F;—8nF, = 0, (12) 
of which the appropriate solution is F,; =8,7* where f; is a constant. 
The equation for f, is unaffected by F, and the equation for F, is 

Fe —39P Fo +57? Fe—9nk, = — 16a, B; 7°, (13) 
with solution 
9 
Fy = 40, 85(n—ge)+Bo 0? = — ihm Bo( = #| +B, 7°. (14) 
The equation for f, is now 
SEVP SEA Saf e—Onfy = — 1604 45 7°+ Op + Fefo— Prof ot Fol i— Ff 
= —2a, Bs(q?—hy* +3897") — 16a, a5 47+G, (15) 
and the solution of (15) with a double zero at the origin is not exponen- 
tially large at infinity if 


a 


0 


- 10 10 2 
G20, P,(0?—" +2 )\r—3 +e dyn = 0 


—titin. 


— 
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from (10). Using Jones’s value for the integral in (10) when r = 6 we 
— 4B 24(4)! a, 8, = —4af, 
whence B, = —ag. (16) 

In principle we can now proceed as far as we like with the asymptotic 
expansion. Extra terms containing logé must be added at the 7th and 
8th stages; the corresponding f’s may be determined in terms of a,, while 
a7, ag depend on a, and a, which now emerges as a new arbitrary constant. 
At the 10th stage however new forcing terms containing (log £)? appear, 
and to cancel these new terms containing (log é)* must be introduced at 
the 9th stage of the asymptotic expansion of y. As a result of satisfying 
the three integral conditions associated with the 10th stage, all the func- 
tions which have appeared up to this stage are determinate in terms of 
three arbitrary constants a,, a;, and a. The development may be con- 
tinued indefinitely and finally the asymptotic expansion of % appears as 
an infinite double series in powers of € and logé whose coefficients are 
functions of » and contain an infinite number of arbitrary constants. 

One interesting feature of these additional terms is that wu will have 


a logarithmic singularity at x = 0. Thus when z is small the additional 
terms in (11) are dominated by 


Bs 


1=s y*logz, (17) 


which tends to infinity as r+ 0. It does not necessarily follow that 
because certain terms, in the expansion of w assumed in this paper, tend 
to infinity as x + 0 for y ~ 0, wu itself is infinite on the line 7 = 0. As 
a counter-example consider Z = y*x”. Then as > 0+, Z—0 for all 
y > 0. On the other hand, for any z + 0 we may expand Z as a series of 
powers of y, viz. 


~ > eet, 
n=0 


every term of which tends to infinity as x > 0+. It would be exceedingly 
difficult, if not impossible, to decide whether u has a logarithmic singu- 
larity at « = 0 by numerical integration from some initial profile. For 
example, Leigh (5) has given a numerical solution correct to six places 
of decimals for a certain boundary layer to within, in his notation, a 
distance 0-00014 of separation. For his values of x = 0-00014 and of y = 2, 
which is just outside the region in which a comparison with Goldstein’s 
asymptotic expansion is possible, the value of (17) above is only 10-5, 
Further, to increase (17) by a factor of ten at the same value of y, we 
should have to be within 10-* of z = 0. The presence of the logarithmic 
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terms makes it clear why Goldstein’s original method of using the separa- 
tion profile as one boundary condition has to be abandoned, for the 
singularities to which they lead make nonsense of the assumption implicit 
in (3). 

The present modification to Goldstein’s theory enables us to resolve 
the difficulties mentioned earlier in this paper. First, it does not matter 
whether the integral condition (10) with r = 6 is satisfied or not. If not 
then logarithmic terms must be introduced at the 5th and subsequent 
stages. If it is satisfied then their appearance is deferred to the 9th stage. 
Second, there is always an infinite sequence of arbitrary constants in the 
asymptotic expansion, which are a,,,, if a, #0 and ay,,, if a, = 0. 
Finally it has been shown (4) that in the related problem of heat conduc- 
tion in one dimension a similar set of arbitrary constants appear in a 
corresponding asymptotic expansion which can be expressed in terms 
of the initial temperature profile. Here therefore it is expected that the 
unknown a’s depend in some way on the initial profile k(y) of u. 


4. The solution downstream of separation (x < 0) 

Although it is not clear whether the extra logarithmic terms make u 
singular on the line x = 0, they render difficult any extension of the solu- 
tion downstream of it. In particular the method suggested by the author 
(6) for a liquid which breaks away from the wall at x = 0 and is there- 
after bounded on one side by a free streamline, cannot directly handle 
these new terms. If the mainstream does not break awayt from the wall 
at x = 0, which it need not do, the only way in which the boundary layer 
could continue downstream of x = 0 seems to be if a, = 0. This implies 
that the mainstream must satisfy some condition at zx = 0 and cannot 
therefore be specified completely independently of the boundary layer. 

If «, = 0, a formal solution for % valid both upstream and downstream 
of x = 0 may be obtained as a double power series in x and y, as Gold- 
stein (1) has pointed out. This cannot be the whole story, however, since 
it implies that the solution in x < 0 is completely specified by U(x), the 
mainstream velocity, and by the solution in x > 0. However, from the 


+ The identification in the literature of the point of separation z, with the point z, at 
which the skin friction vanishes may lead to some confusion. In this paper 2, is the earliest 
point at which the fluid moving forward in the boundary layer is no longer adjacent to the 
wall, but is separated from it by a layer of fluid moving in the reverse direction. It coincides 
therefore with x). However, the word separation is sometimes used to describe the point 
£, at which the main stream is observed to leave the neighbourhood of the wall. In this 
paper to distinguish the two phenomena we use the term breakaway to describe zp. There 
is no satisfactory theory at present which can predict where breakaway occurs. It may 
be observed to occur even if the skin friction does not vanish and conversely the mainstream 
has remained near to the wall even when the boundary layer has separated. If breakaway 
and separation both occur then z, is usually just downstream of x, = 25. 
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form of the boundary layer equations it follows that disturbances are 
propagated by convection, with velocity u, in the xz direction and by 
diffusion in the y direction. Thus upstream of separation (x;,y;) can 
only be modified by a disturbance at 2, y if x >2z,. Downstream of 
separation the boundary layer can be divided into two parts in one of 
which w > 0 and in the other u <0. Hence wherever in x < 0 the 
disturbance is introduced it diffuses into both parts of the boundary layer 
and then spreads by convection into the whole boundary layer downstream 
of separation. The boundary layer upstream of separation is of course 
unaffected except indirectly through any consequent modification of the 
mainstream. 

Mathematically this means that in 2 < 0 extra terms, additional to the 
double power series, must appear in the solution, which are not fully 
determinate by the conditions in the neighbourhood of z = 0. A method 
is given here by which these terms may be found, and it will appear that 
at x = 0 not only do they vanish but so also do all their derivatives with 
respect to x; thus the join between the flows upstream and downstream 
of separation is perfectly smooth. 

Let the stream function obtained by the expansion in a double power 


series be ys, equal to hy?+- 4a, ry?-+Jo (18) 


near x = y = 0 where a, is a positive constant and g, is Ofy*(x*+-y*)}. 
Further let the complete stream function » = %)+y, where ¢, <p 
when x is small. Near x = 0, therefore, 4? may be neglected in com- 
parison with %,, which, on substituting into (1), satisfies 


2, Oly My _ OMby Mths , Hy Pho Ar Mo _ gy ag) 


oy ba Oy Oxdy Oy? Ox ~ Oy Oxdy dy® Ox 


with boundary conditions 
uv, — “1 — 0 when y = 0,2 <0, %, = Owhenz = 0,y > 0. 
cy 


We look for a solution of the form 

%, = ‘Vi(0,x), where 0 = —y/a,2, (20) 
assuming that while é‘¥,/éx may be much larger than ‘f,, é'f,/20 is of the 
same order as ‘Y,. This assumption may be verified a posteriori because 
‘Y, has an essential singularity at x = 0. The equation satisfied by , is 
ay, ov, 
Aap %3 109-1) + Oz) = O(z' ¥}) 

(21) 


where the terms O(z) are of the form x multiplied by a double power series 


A + dag 200-2) + 0(2)] 
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in x, @ and the terms O(2°‘V,) are of the form z* multiplied by a double 
power series in x, 6 and either ‘YY, é'Y,/26, or é,/26. The simplest 
appropriate form for ‘Y,, but not oe the only one, is 


¥,= exp: B ro ay > arn) (22) 


; | 3a, a9 cai 2a x? =| 


where f, y, 5, n are constants of which 8 must be positive and F,(@) is 
a function of @ only satisfying the boundary conditions 


F(0) = FX(0)= 0 and F(@)++ 0 exponentially as 8-0. (23) 
The last condition is necessary because Y, must be bounded as y > © and 
because the disturbances which affect ‘4, are only propagated in the 
direction of x increasing near the plate and reach large y only through 
diffusion. The equation satisfied by J, is 
Fy — $80(0—2) Fo +B(O—1)Fy = 0, (24) 
which is too complicated to lead to a simple solution. Here therefore 
it will only be shown that solutions exist for an infinite number of 
acceptable values of 8 and that then (22) contains three arbitrary con- 
stants. Differentiating (24) with sey to 6 il 
Fi —3B6(0— )Fo+BR, = os 
We now suppose that f is large so that BF, may be neglected in com- 
parison with Fly. Thus the problem is reduced to solving 


iv_489(0—2)F* = 0 


subject to the boundary conditions F4(0) = 0, Fj ++ 0 exponentially as 
6» co. The range of @ has to be divided into four parts and each con- 
sidered separately: 


(a) @ > 2 and f'(@—2) large. Here the appropriate solution is 
yi 
Fy ~ [(6—2)]* exp) — Bt | [6(8—2)}! ad 
2 


where ~ means that the ratio of the left- and right-hand sides tends to 
a constant value as 8 + oo. Then near @ = 2 
", ~ (0—2)-4exp{—(28)(0—2)}}. 
(b) p'(@—2) = O(1). Here to match (25) 
Fo ~ BHO—2)'K {3 (2B)'(8 
and when f#(@—2) is large and negative 


Fy ~ (2—86)-* sin{ fa — §(28)#(2—6)!} 
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(c) 0< 6 < 2; BO, (2—6)Bt large and positive. Here, to match (26), 


6 
Fw [0(2—6)]- sin} [ (0(2—0)}§ 40—4nB!+ 4 


and near 6 = 0 Fo ~ 9-* sin{}a—4aBt+ §(28)108}. (27) 
(d) 6st = O(1). Here since Fj (0) = 0 
Fo ~ BOY {3 (28)'0}, 


and when 68? is large 


| 0-sin|§(2p)40 -= : (28) 
Hence, to match (27) and (28), 
en — =" = }n—4nfi, 


where s is a large positive integer. Thus 
B = (3s—1)? (29) 
and there are an infinite number of suitable values of f. 

Each such 8 generates an independent solution of (21), of the form (22), 
containing three arbitrary constants, which may be taken to be F%4(0), 
F‘(0), F3(0). Given these in addition to £ all other terms in (22) may be 
found. First y is found from the equation for F,(0), which is 


FY —}80(0—2)F,+B(0—1)F, = yFo+G,(0), (39) 


where G, is a linear function of F, but is independent of y. Since the 
solution of the homogeneous equation for F, which satisfies the conditions 
at @ = 0 also satisfies the condition at 6 = oo, a solution of (30) satisfying 
these conditions is only possible when 


@ 


i F?(yF” +-G,) dd = 0, (31) 


0 
which determines y in terms of 8 only. Since an arbitrary multiple of F, 
may be added to #,, F}(0) is arbitrary. From the equation for F,, 8 may 
be found in terms of F5(0), Fj(0), while F3(0) is arbitrary. Similarly n 
is determined from the equation for F, and F7(0) (r > 3) from the equation 
for F,,,, thus formally completing the solution. 

The form of ‘Y, obtained above is very acceptable, for it has an essential 
singularity at x = 0 which ensures that it, and all its derivatives with 
respect to x, tend to 0 as x > 0—, so that the join with the solution in 
x > 0 is perfectly smooth. Further, as would be expected, because there 
are an infinite number of suitable ‘, containing in all an infinite number 
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of arbitrary constants, the solution in x < 0 is capable of an extremely 
wide variation depending on the boundary conditions beyond separation. 
The situation is similar to that in the theory of the impulsive motion 
of a flat plate in an incompressible fluid studied by the author (7). There, 
if x denotes distance from the leading edge of the plate, ¢ time from the 
start of the motion, and U the constant velocity of the plate, a similar 
difficulty arises in the region 0 < € < 1 where é = 2/Ut. The stream 
function ys is known at € = 0 and at & = 1 but a step by step integration 
cannot be started from either point because the flow depends on boundary 
conditions at € = 0 and at = 1. Essential singularities develop at 
€ = 0, 1 containing an infinite number of unknown constants which 
depend in some way on the conditions at € = 0, 1. There is a close 
similarity between the essential singularities at € = 0 and n the present 
problem, the forms being alike and both being centred on y = 0. 
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SUMMARY 
Similarity solutions describing the flow of a perfect gas behind strong shock waves 
are investigated for the three cases of plane, cylindrical, and spherical symmetry. 
The flow is caused by an expanding piston, and the total energy increases as a power 
of the time. The ratio of kinetic to internal heat energy of the gas is computed and 
it is found that in the case of a piston which is expanding at a uniform rate there 
is equipartition of energy. 


1. Introduction 

In this paper a class of exact solutions of the equations governing the 
motion of a perfect gas is given and the properties of the solutions 
discussed. The flows depend on one spatial coordinate only. While certain 
of these solutions have been studied at different times (1, 2, 3, 4), it is 
nevertheless of interest to examine a general set of such solutions in order 
to gain a better understanding of similarity flows; in particular the problem 
of the occurrence of singularities in the flow and the fitting of appropriate 
boundary conditions is investigated in some detail. The solutions to be 
discussed are of the progressing wave type, or similarity solutions. These 
are obtained by assuming the solutions to be of a certain form, which is 
given below, and by virtue of this assumption reducing the mathematical 
problem to one involving ordinary, instead of partial differential equa- 
tions. This simplification is achieved by postulating that the dependent 
variables in the problem—the pressure, density, and velocity—can each 
be written as the product of a function of the time, and a function of 
a single variable x which is of the type 


xz = rt-, (1.1) 


where r is the distance from the origin, ¢t is the time, and A is a constant 
which is determined, effectively, from considerations of energy. Such a 
set of solutions can then be employed to describe the flow behind an 
infinitely strong shock which is advancing into a gas at rest. It is possible 
moreover, (4), to treat the case in which the density distribution ahead 
of the shock front has the form 


p= per, (1.2) 
+ Present address: Royal*Military College of Science, Shrivenham. 
[Quart. Journ. Mech. and Applied Math., Vol. XI, Pt. 4, 1958] 
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where « is a constant. When the density of the undisturbed gas is uniform, 
and it is assumed that the total energy carried by the flow remains 
constant, the resulting solution, first obtained by Sir Geoffrey Taylor (2), 
describes the initial stages of an intense explosion. Now the total energy 
in the flow is the sum of the kinetic and internal heat energies of the gas; 
in practice there will be losses due to dissipative effects while there will 
also be a gain in the internal heat energy as the shock front advances and 
encloses more of the quiescent gas. However, this increase has effectively 
been omitted because of the assumption that the shock is infinitely strong, 
which is equivalent to assuming that the pressure, and therefore the in- 
ternal heat energy of the gas ahead of the shock front, is negligible. In 
the present paper it is proposed to investigate the similarity flows which 
exists behind infinitely strong shock waves when the restriction of constant 
energy is removed, and the total energy is allowed to vary with time in 
the manner E = Ext’, (1.3) 
where £ is the total energy and £, and s are constants. This class of 
flows includes the blast wave, which is given by the value s = 0. Further- 
more, attention will be confined to positive values of s only, that is, to 
those cases in which the total energy increases with time. Since the motion 
of the gas is adiabatic, and the shock is infinitely strong, this increase can 
only be achieved by the pressure exerted on the gas by an expanding 
surface—plane, cylinder, or sphere. Thus the flow is headed by a shock 
front and has an expanding surface as an inner boundary. The position 
of this inner surface will be determined from integration of the equations 
of motion; accordingly a step-by-step method of integration is used, 
starting at the shock front and continuing into the gas until the surface 
is reached. Particular attention is paid to the energy of the flow: in the 
blast wave nearly 80 per cent. of the energy present is in the form of 
internal heat energy and it will be seen that one of the effects of an in- 
crease in the value of the parameter s is to reduce this ratio, until in the 
case of a uniformly expanding surface there is an equipartition of energy. 
This result is of some interest in view of recent applications of similarity 
solutions to problems in hypersonic flow (5). Throughout the subsequent 
working it will be assumed that the density of the undisturbed gas ahead 
of the advancing shock is uniform. Numerical solutions for spherically 
symmetric flows have been obtained by Sir Geoffrey Taylor in the two 
cases s = 0 (2) and s = 3 (1). The solution for cylindrically symmetric 
flow with s = 0 has been found numerically by Lin (6), while analytic 
solutions with s = 0 in the three cases of plane, cylindrical, and spherical 
flow have been noted by Sakurai (7). 
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2. Equations governing the motion 


The equations describing the one-dimensional unsteady flow of a perfect 
gas are 


(2.1) 
= FON. FOP PN 
at" optP pte =o (2.2) 


< (pp?) +u (pp) = 0, (2.3) 


where wu, p, p are the velocity, pressure, and density of the gas, and k takes 
the values 0, 1, 2 for the respective cases of plane, cylindrical, and 
spherical symmetry. In the subsequent analysis it will be convenient to 
replace equation (2.3) by the equivalent equation 

Of, 21? \,) Af 2, YP \\_9 24 

= (do +P )+5 = lag }pu +r * (2.4) 
Equations (2.1), (2.2), and (2.4) are the fundamental set to be integrated, 
subject to boundary conditions which will be discussed below. Turning 


now to the choice of a similarity variable, it will be recalled that in the 
case of the blast wave problem this variable was 


z=r/R, (2.5) 
where RF is a function of time only and gives the position of the shock 
front at any instant, and this definition will also be employed in the 
present work. Accordingly, the velocity of expansion of the shock front is 

V = dRidt (2.6) 
and the shock front is represented by x = 1. Solutions of equations (2.1), 
(2.2), and (2.4) are now sought in which the velocity, density, and pressure 


have the forms u = Vf(z), (2.7) 


Papel em (2), (2.8) 


p = poh(z), (2.9) 


where p, is the density of the uridisturbed gas ahead of the shock front. 
The total energy of the flow may be written as 


iu { Jpu? dr+ | P ; dr, (2.10) 
y— 

where dr is an element of volume; the first integral represents the total 
kinetic energy of the gas, and the second integral the internal heat energy. 
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2. Equations governing the motion 
The equations describing the one-dimensional unsteady flow of a perfect 
gas are 
oO ig Oe ce (2.1) 
a or por’ 
pa = 


ep ep at 


< (pp) +u rile = 0, (2.3) 


where u, p, p are the velocity, pressure, and density of the gas, and k takes 
the values 0, 1, 2 for the respective cases of plane, cylindrical, and 
spherical symmetry. In the subsequent analysis it will be convenient to 
replace equation (2.3) by the equivalent equation 

i | ee 1 ele 2, YP \\ _ 9. 2.4 

(tou +P) +5 S[r*u( dour + 22) — (2.4) 
Equations (2.1), (2.2), and (2.4) are the Een set to be integrated, 
subject to boundary conditions which will be discussed below. Turning 
now to the choice of a similarity variable, it will be recalled that in the 
case of the blast wave problem this variable was 


x=r/R, (2.5) 
where F is a function of time only and gives the position of the shock 
front at any instant, and this definition will also be employed in the 
present work. Accordingly, the velocity of expansion of the shock front is 

V = dR/dt (2.6) 
and the shock front is represented by x = 1. Solutions of equations (2.1), 
2.2), and (2.4) are now sought in which the velocity, density, and pressure 


have the forms u = Vf(z), (2.7) 


p = Po" (2), (2.8) 
Y 


p = poh(z), (2.9) 


where p, is the density of the undisturbed gas ahead of the shock front. 
The total energy of the flow may be written as 


(2.10) 


where dr is an element of volume; the first integral represents the total 
kinetic energy of the gas, and the second integral the internal heat energy. 
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By using the assumed forms for the dependent variables this can be put 
in the form 


1 
E = pyc, V2RE [ (yng24—2 or ae, 
J \ yv(y—1)) 


where «€,, = 2*z#*3-”, and 2, is the coordinate of the expanding surface. 
Since f, g, and h are functions of x only, this integral will be a function of 
k and y, and (implicitly) s. Now the type of flow under consideration is 
one in which the energy increases with time, and the law of variation has 
already been stated, namely equation (1.3). From this relation and 
equation (2.11) it follows that the motion of the shock front is described 


by the equation 
: ; dk pen | Ey |e 


Rue 1) 
dt €, Pod 


To 


where Fx { f nf? +555 |e dx 


aa) 


Equation (2.12), on integration, yields 


* (: eo = a (Seg) ene, (2.14) 


+2 €x Pod 


where it is assumed that R = 0 initially. From the relation (2.14) it can 
be seen that the value s = k+-1 corresponds to the uniform expansion 
of a plane, cylinder, or sphere; the solutions of physical significance appear 
to be those for which s lies in the range 0 to k+1. The basic set of equa- 
tions (2.1), (2.2), and (2.4) can now be put in the form 


s—k—1 


(e—f)f' =: iS ee we aan 


airs 


y—1)} 


1- 
2(k+- -8) 2* EB, (2) 
8+2 
where EB, = jf? —_2_ 

vy—1)’ 

and a prime denotes differentiation with respect to z. The boundary 
conditions must now be examined and these are discussed by Kynch (3); 
in the present notation, for very strong shocks, these are 


9 
shaders (2.18) 
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2y 
1) = —“, 2.19 
g(1) mee ( ) 


MU aot. (2.20) 
y—1 
Equations (2.15), (2.16), and (2.17) are now integrated numerically using 
the Runge—Kutta-Gill method starting at the shock front and continuing 
until a value x, is reached such that f(z,) = 2). This is simply the 
kinematic condition at the expanding surface which states that the velocity 
of the gas is equal to the velocity of the surface itself. In view of the 
similarity assumption concerning the flow, the total mass of gas between 
the expanding surface and the shock front should be the same as that 
originally contained in the volume given by 0 < x < | and this fact was 
used as a check on the accuracy of the numerical integrations. The results 
of the integrations are shown in the following tables; the values of k and y 
used are indicated for each table. The second column gives the value 2, 
at which the surface occurs, and the third column the value of the pressure 
at this point. The fourth column gives the value of the integral J defined 
in equation (2.13); this is used in evaluating the total energy carried by the 
flow. The integrations were performed on the University of Illinois digital 
computer. j 
TABLE 1 


y=14. k= 2. Spherical flow 





Ss X» (= f(x) (9) J 


° ° 04264 0°4264 
0578 04854 0° 3960 
o2 0-736 05748 0°3610 
os 0846 0°7302 o°3173 
1°00 0"g00 00-9079 0°2832 
2°00 0-932 Ir215 02556 
0-942 1*2450 02526 

















TABLE 2 
k=1. Cylindrical flow 





me | Bae J 


° 0°4351 o-6414 
0°408 o-4814 0°5900 
0468 04986 0°5773 
05232 0° 5607 
0°5975 05178 
07206 04660 
o-8211 0°4344 
09776 03992 
1°2246 0°3702 
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TABLE 3 
y=14. k=0. Plane flow 





S Xo &(Xq) 





° ° 0°4550 
0-221 0° 5086 
0365 "5591 
0°536 0°6523 
0°694 0°8139 
0°736 0°8846 
0°768 "9500 
0808 1°0662 
0°833 1-1667 














TABLE 4 
= 2. Spherical flow 





g(%9) J 


04621 08592 
0°657 o°5055 07068 
o8r1 05749 0°5637 
0°903 06999 0°4309 
0942 08458 0°3541 
0-962 T°O251 0°2949 
"970 1*1275 0°2882 

















TABLE 5 
k=1. Cylindrical flow 





Xe (Xo) J 





00-4648 1°2892 
501 05006 1'0785 
567 O°5142 1°0203 
637 0°5338 0-949! 
762 0°5938 0°7932 
857 06951 0°6407 
8096 0°7786 05644 
929 O-gIOI 0°4860 
"954 1°r182 o°4251 


2209999900 














TABLE 6 
k=0. Plane flow 


&(X9) J 

O-47I5 2°5200 
05184 19722 
0° 5626 1°6699 
06440 1°3439 
07849 1°0653 
08464 09959 
0°9033 0°9365 
10043 08640 
T0909 o-8271 
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TABLE 7 
Values of go for different values of y’ 

















y Plane Cylindrical | Spherical 
ro 050000 0*50000 0" 50000 
vl 0°48361 048158 0° 48083 
1°2 O°47154 0°46476 046211 
I°3 046231 0°44932 0°44401 
1°4 0°45502 0°43507 0°42637 
1°67 0°44149 o-40188 038274 
20 0°43124 0° 36788 0°33442 
30 0°41693 0°29630 0°20233 





3. Discussion of results 

Although the integration of equations (2.15), (2.16), and (2.17) was 
performed numerically, it is of interest to note that a partial analytic 
solution of the problem exists. Writing equations (2.3) and (2.16) as 


TW ip DZ), 8) 
bp +2) 4, (3.1) 
cee’ 
c~s' ts} (3.2) 
a first integral can be obtained by multiplying (3.2) by 
2(s—k—1) 
y+ REED 
and subtracting the result from (3.1); the resulting equation can then be 


integrated, yielding g = C{xk(x—f)} hr, (3.3) 


where C is a constant which can be determined from the shock front 
boundary conditions, and 


_ 2(1+k—a) 
~ (8-+2)(k+1)" 
Inspection of equation (2.17) shows that this equation integrates imme- 
diately if 2(k-+-1—s) . 
= k+1, 
s+2 
or (k+3)s = 0. (3.4) 


This leads to the analytic solution for the blast wave in a uniform atmo- 
sphere. A complete analytic solution also exists in the case when s = 1 
and k = 0, which represents the flow in front of a plane piston moving 


with a uniform velocity into a gas which is at rest. The expressions for 
5092.44 Ee 
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the functions f, g, and h are 
9 


2) = —, 3.5 
fle) = (3.5) 
9 

et 
ae. ag 3.6 
g(x) ae (3.6) 
hx) = YEE, (3.7) 

y—1 
The value of x, which gives the position of the piston relative to the shock 
front, is %y = 2/(y+1) (3.8) 
and the expressions for the velocity, density, and pressure of the gas are 

oV 
‘= 27" ’ (3.9) 
y+1 

l 
p= port, (3.10) 

y—l1 

2p, V? 
i 3.11 
Pa (3.11) 


where V is the (constant) velocity of the shock front; the piston is advancing 
with a speed of 2V/(y+1),so that the volume occupied by the flow steadily 
increases with time. The total energy of the flow contained in a volume 
having unit cross-sectional area is 


1 
é=pveR | [ape -| d 
y 


a1 (y—1) 
2/(y+1) 
1 
, 4 4 4p,V?R 
== 2 ae d: => ae aT © vo. 2 
reves) | * +I? ssi 


2/(y+1) 

It will be noticed that the values of these two integrals are the same so 
that the kinetic energy is equal to the internal heat energy of the gas in 
this case. Furthermore, the rate at which the total energy of this volume 
is increasing is 3 

dé _ 4p) V* (3.13) 

dt (y+1)? 
Now this increase is equal to the rate # at which work is done on the gas 
by a unit area of the expanding piston; this latter quantity is simply the 
product of the surface pressure and the velocity of expansion of the plane, 
so that ie 2po de 2Ve ae 4p, v3 (3 14) 

eth yt (y+) 

Returning to the general case in which a complete analytic solution does 
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not appear to exist, the rate at which energy is being fed into the flow 
is again equal to the product of the surface p: assure and the area and the 
velocity of the expanding surface so that 


R = (Xp Ry Pog ayV (3.15) 


ae 


Cylinder 








1 ¥: 
Fic. 1. The variation of the integral J, defined by equation (2.13), with the index a. 
In each of the three cases, only the values of s in the range (0, k+-1) were employed. 
The value of y used in these calculations was 1-4. 





Now the total energy of the flow is given by equation (1.3) and hence 
R = sk, t, (3.16) 
If these equations are now combined with (2.14) it follows that the integral 
J can be expressed as ‘ 
2) 

J = ET?) perigye 3.17 
provided that s does not vanish; it is found that there is good agreement 
between the values of J determined from the numerical integration and 
those obtained from (3.17), using the computed value of the surface 
pressure of the gas. 

The manner in which J varies with s is shown in Fig. 1, where the three 
cases—plane, cylindrical, and spherical—have been computed using a 
value of 1-4 for y. The figure shows that in each case the value of this 
integral decreases as 8 increases from zero to k+1. At first sight this 
appears to be somewhat surprising, and it is of interest to investigate 
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this integral a little more closely. Now the energy is comprised of two 
parts, the kinetic energy and the internal heat energy of the gas. When 
s = 0 the contribution of the heat energy is far greater than that of the 
kinetic energy; in the spherical case, Sir Geoffrey Taylor (2) has calculated 
these quantities for several values of y. When y = 1-4 the kinetic energy 


% KE 
4 











o ARS 6 enema. «et, 
Fic. 2. The variation with s of the kinetic energy of the flow, expressed as a per- 
centage of the total energy. The spherical case only is shown for the values of y 
indicated. Note that as s approaches the limiting value of 3, the fraction approaches 
50 per cent. so that there is equipartition of the energy. A similar result holds for 
the plane and cylindrical cases, the limiting values of s being 1 and 2 respectively. 


represents only 21-8 per cent. of the total energy present and for smaller 
values of y this fraction is even less. This is due to the fact that the 
pressure tends to a constant value as x + 0, whereas both the density and 
velocity tend to zero at the centre so that there is very little contribution 
to the kinetic energy in the region 0 < x < 0-6. In the present case, the 
region in which the flow is occurring is smaller than that in the blast wave 
and since the pressure does not rise to an extremely high value compared 
with its value immediately behind the shock front, the contribution of 
the internal heat energy will be considerably less than in the constant 
energy solution. On the other hand the contribution from the kinetic 
energy integral rises with increasing values of s; but this increase is 
not sufficient to make up for the loss in the heat energy integral. Fig. 2 
shows the variation of the kinetic energy expressed as a percentage of the 
total energy, again for spherically symmetric flows only. It will be seen 
that there is a sharp rise in this fraction as s increases from zero through 
small positive values, and then the ratio increases more slowly up to 
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s = 3 when there is approximately an equipartition of energy. The com- 
puted value is slightly over 50 per cent., but this is probably due to 
accumulated errors of integration. In the case of plane flow with s = k+-1 
the analytic solution displayed above shows that the two energy integrals 
are equal and a similar result may be expected in the spherical and 
cylindrical cases. It is worth noting at this point that in the well known 
solution of Primakoff (8), which describes the motion of a spherically sym- 
metric blast wave in water, there is also an equipartition of energy. 
Furthermore, in the case of blast waves advancing into non-uniform 
atmospheres (4) the equipartition of energy is associated with the onset 
of cavitation in the flow: when no cavitation is present the kinetic energy 
is less than the internal heat energy in the flow, but when cavitation occurs 
the reverse is found to be true. For the critical flow, which separates these 
two regimes, there is an equipartition of energy. 

In the constant energy solution, the blast wave, the pressure of the gas 
at the centre can be computed from the analytic solution obtained from 
equations (3.3) and (2.17) where s is put equal to zero. After some reduc- 
tion it is found that the expression for the central pressure can be written 
as 


Jo = YR DARD A-D)Q AIR) (yy 4 J R-II R+D-D) > 


—(k—1)y+3k+1)0-® 
2y+k—1 é 
ipa (k?+-2k+5)y?+ (—3k?+ 2k+ 1)y+4(k?— 1) 
(k+3)yk+1)—(k—1)} 

The values of g) can then readily be calculated and the results are shown 
in Table 7; some of these values are also given in Tables 1-6. The results 
show some points of interest: in the first place, it will be seen that in the 
case of isothermal flow the three values of g, are equal. As the value of 
y increases, the central pressure gradually falls, and as would be expected 
the decrease is greatest for the spherically symmetric flow and least for 
the plane flow. 

The density of the gas at the expanding surface is next examined and 
the partial analytic solution (3.3) is of use in this connexion. In the blast 
wave the density vanishes at the centre and the density gradient is also 
zero there. For the present case, equation (3.3) shows that, in view of the 
fact that the pressure remains finite at the surface, 


h(x) ~ (x—f Py, (3.18) 


When s = 0 this index has the value 1/(y—1) and when s = k+1 it 
vanishes; for intermediate values of s it is positive. Consequently, except 





where 
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for the value s = k+-1, the surface density of the gas is zero. Moreover, 
since the pressure is finite, the surface temperature will be infinite while 
the density gradient will be zero, finite, or infinite depending on whether 
) 
y—5 
or, expressing this in terms of s, 
a 2(k+1)(2—y)_ 
© 4-+y(k+1) 
When s exceeds this critical value, there is accordingly a very narrow layer 
on the surface across which the density falls abruptly to zero, the tem- 
perature rises to infinity while the pressure remains finite. In practice, of 
course, the effects of viscosity and heat conduction are of importance in 
this boundary layer and the similarity solution would break down at some 
point in this region. 


$1 


(3.19) 
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SUMMARY 


In the present paper, the transverse component of the velocity in the laminar 
motion of a plane symmetrical jet of a compressible fluid has been found, and its 
variation with respect to the distance from the jet axis has been studied. 


1. Introduction 

D. G. Tooss, in a recent paper (1), has obtained a complete closed 
analytical solution for a plane compressible symmetrical jet under the 
usual boundary layer conditions. From his expression for the transverse 
component of the velocity, he has deduced that this transverse velocity 
is always directed towards the jet axis. He has also obtained an equation 
to determine the stream function and from this we can also obtain the 
transverse component of the velocity. We find, however, that this ex- 
pression is different from his and, on a closer examination, that his first 
derivation of the expression for the transverse velocity is incorrect. We 
have found the correct expression for this and it appears that, for any 
given section of the jet, the transverse velocity is directed away from the 
axis of the jet at points near the axis. As the distance from the axis 
increases, this velocity component first increases and then decreases, till 
at some definite distance it becomes zero. After that the transverse com- 
ponent is directed towards the jet axis. The present conclusions appear 
also to be physically more plausible thai. the conclusions of (1). The 
streamlines and the velocity profiles turn out to be of the same type 
as fer incompressible jets studied by Bickley (2). In fact, while in this 
degenerate case, the formula for the axial velocity component wu in (1) re- 
duces to the corresponding formula for wu in (2), the formula for the trans- 
verse velocity component v in (1) does not reduce to the corresponding 
formula for v in (2). We verify that our formula for v reduces to Bickley’s 
formula in the degenerate case. 

We use the notation of (1) throughout. 
(Quart. Journ. Mech. and Applied Math., Vol. XI, Pt. 4, 1958] 
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2. Derivation of transverse component of velocity from the stream 
function 


From (19), (25), and (38) of (1) we obtain 


ii 
y =v I. te) = Wr X#lo (ad 





+Hy—1)ag| Y— —wx4 2% [+4072 (1) 

Also 
6. eS. 2. fee (2) 
u pu — —p,(eb/ey) ss OW Jey V(u;/v; LO)ev /eY” is 


From (1) and (2), 


u;l\v T, v6+¢€ 2 re 
Je oa Hil ys" ye—¢ oa] HHO ge F]. 


But, from (1), 





u = u,X—4*(1—}3l*); (4) 


ee ee u;v;C\(T, 2 $F ss 
hence v= Jt 7 Vials val! 6 log = t Cj}+ 


hg 2 . 
+My MI (1- -4) \ (5) 





On the other hand, from (37) and (40) of (1), we get 
2 
‘ v= — (4 5 Nirx-| 7, —4(y—1) MP X (1-2) + 


+0 DMJ—2}X- —*}} (6) 


Expressions (5) and (6) are obviously different, and we examine the reason 
for this discrepancy and its consequences in the next section. 


3. Earlier derivation of v examined 
Since di = (=). dx +. +(4)_ a x () ay +(“) dx 
ve eG ey Ox} y 
it follows that 


fae /(“% v5 \x- 44 XK (u,v, Lt r(¢ =) | de Xu a 1 | dy 


z 


— OM dy—P? de. (7) 
P; Pj 
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From (7), 0 = —afat, |S) x44 Huy», L0)(2) | (8) 
and ve al Xu % wo) |. (9) 


Instead of (8), (1) gives 


v= — 5] (“HES (10) 


which is obviously incorrect. In (1), (10) has been used tou determine v 
and the boundary conditions for solving the basic differential equation 
and to deduce that the transverse velocity in the jet is directed towards 
the jet axis. We examine the necessary modifications below. 

From (13), (16), (18), (19), and (25) of (1) 


a | is ae, SE +7]: 
(ee), ~ Tee EO 5p PPE sty 
Substituting (11) in (8) gives 
{ = YZ-4. (12) 
Using (19) and (25) of (1) we see that (8) is an identity which cannot 
determine v. The only way of determining v is therefore through the 
stream function, as has been done above. 
In (1), (10) has been used to deduce the boundary condition 
€=0 when y=0. (13) 
Fortunately (13) can also be deduced from symmetry considerations, for 
by symmetry y = 0 is a streamline and since numbering of streamlines is 
arbitrary, we can take 
Y=0 when y= 0. (14) 
(13) then follows from (12) and (14). 
From (9), we also get 


v(v; LC) 
(4) = hun) — MO xr. (15) 
Integrating, ~ eee dl+¢4(X), (16) 


where @ is an arbitrary constant al ¢(X) an arbitrary function of X. 
However, since from (13) we have y = 0 when { = 0, it follows that 


¢(X)=0, a=0 


and so y= Je) Xi j ¥O dé. (17) 
0 


This is the same as the corresponding equation in (1), but it is seen to be 
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true, not only when X is held constant, but also when its variation is 
taken into account. 


4. Discussion of variation of transverse component 

In (5), the coefficient of 1/X* is always positive, while the coefficient 
of 1/X* is positive for small values of f and takes negative values for 
larger values of ¢. Also v = 0 when (¢ = 0 and wv approaches a finite 
negative value when ( > v6. Again, from (15), (@y/@C), is always positive, 
so that for every fixed x, as y increases { increases. Thus, for every x 
there is a value of y such that up to that value, the transverse velocity 
is directed away from the jet axis and after that the transverse velocity 
is directed towards the jet axis. 

The equation of the curve which divides the regions in which transverse 
velocity is directed towards or away from the jet axis is given by 


| (4p) = Text toe(F*2) —1—nag(s—F) +x, —72), 
I \v, 1 6—f 
(18) 


T l v2 5+ F l c2\2 
o = <2|1(; en, y—1)M? _o\" 9 
xi [a v6- came 3X1” ( i) a 


where ¢ is the parameter. 
For a fixed X, putting (év/éf) = 0, we get the value of ¢ for which v 


has an extreme value and then (1) and (11) determine the catalan 
value of y. 


5. A degenerate case 
In the degenerate case, when ¢; = t,, and the orifice velocity u,; is very 
small, proceeding as in (1) we find 


y 


as = tanhé, 
v6 
where £ — XY i OM $ 
; 3x 16p; v} 
In this case (5) gives 


X ca. 7) - X-tsech2£, 
Se 


ee V3 X 
which agrees with Bickley’s formulae. 


5% 


i{ {2é sech*é—tanh €}, 
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SUMMARY 


Variational principles for the approximate computation of the elements of the 
scattering matrix for anisotropic obstacles in waveguides are presented. 


1. Introduction 


THE effect of a ferrite obstacle in a waveguide on the propagating modes 
of the guide may be described in terms of the scattering matrix, which 
relates the reflected and incident amplitudes of the propagating modes 
at suitably chosen reference planes. The calculation of the elements of 
the scattering matrix involves a knowledge of the electromagnetic field 
within the obstacle (equations 16-21), which has to be obtained from a 
solution of Maxwell’s equations everywhere within the guide subject to 
appropriate boundary conditions. 

The exact solution of this problem is a formidable task if not an im- 
possible one. Quite naturally therefore, we seek approximate solutions. 
This is the purpose of this paper. 

The method which we develop here is an extension of the work of 
J. Schwinger (1) on isotropic obstacles in waveguides to anisotropic 
obstacles. With the introduction of an appropriate dyadic Green’s func- 
tion we are able to obtain a formal solution to the problem in terms of 
several integrals involving the field vectors within the obstacle. The 
resulting integral equations are by no means easier to solve. It is possible, 
however, to set up variational expressions for the elements of the 
scattering matrix, which are stationary with respect to first-order varia- 
tions of the fields about their true value. 

Consequently, we have a very powerful method for obtaining approxi- 
mate solutions for the elements of the scattering matrix. 

+ The research reported in this document was supported jointly by the Army, Navy, 


and Air Force under contract with Massachusetts Institute of Technology. 
t Staff Member, Lincoln Laboratory, Massachusetts Institute of Technology. 


[Quart. Journ. Mech. and Applied Math., Vol. XI, Pt. 4, 1958] 
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2. Integral equations for the electromagnetic field and elements 
of the scattering matrix 
We consider the problem of an obstacle with tensor electromagnetic 
properties within a waveguide. The electromagnetic fields within the 
guide, assuming a time factor e’”, satisfy the differential equations 
VxE= —jopoH—jop'.Hy (1) 
and VxH = jwe, E+jwe’.E 
€,) and py are the electric permittivity and magnetic permeability of the 
empty guide. The tensor magnetic permeability p and the tensor electric 
permittivity € of the obstacle are given by 


B=’ +pol, e= e'+e,I, 
where I is the unit dyadic; e’ and pw’ in equation (1) are to be taken zero 
outside the obstacle. 
The wave equations satisfied by the electromagnetic field vectors are 


VxVxE-—FE = wp, €’. E—jwV x (p’.H) (2a) 
and VxVxH—F?H = we, p’.H+jwV x (e’.E), (2b) 
where k? = we, pio. 


To solve equation (2a) we utilize the electric dyadic Green’s function 

N(e,z\p’,2’), which satisfies the differential equation 
VxVxN—FN = [8(r—r’) (3) 
and the boundary condition 
nx N(p,z!e’,2’) = 0 

when e = xi+yj lies on the boundary of the guide. At the reference 
planes N must represent radiated waves moving in the increasing 
z—z'| direction. If the electric field within the obstacle is non-divergent 
then we shall, in addition, require the electric dyadic Green’s function to 
be non-divergent too. Such is the case, for example, for ferrite obstacles 
for which e€ = el. 

The non-divergent part of the dyadic Green’s function N, is expansible 
in terms of the electric fields of the usual TE and TM modes of the guide. 
If E{Ye-%y™ is the electric field of the nth mode moving in the increasing 


z-direction, and E{e/y"* the electric field of the wave moving in the 
decreasing z-direction, we find that 


a E‘)(9)E2)(9’)e-ivme-2) (2 > 
N,(e,2/ 9’, 2’) _ >) A, x n (pe) n (e’) : ( - 
n= ED (p)ED(e' jem (z <2 
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where 


~f ( Ke 2 (/w "ug [ H? dS (for TE modes) 
" 2yn | k&, /k* { H2dS (for TM modes) 


The curl-less part of the dyadic Green’s function N,, which we should 
expect to be derivable from the gradient of a vector function, satisfies the 
differential equation 


(ken = Ky). (5) 


—kY.N, = Vi(r—r’). 
It is readily shown (2) that 


N,(rir’) = — ,VV'G(r\r’), (6) 


where G(rir’) is the scalar Green’s function satisfying the differential 
equation VG = —a(r—r’). (7) 
In order that nN, should satisfy the boundary conditions, G(r|r’) 
should be zero on the boundary of the guide and likewise be exponentially 


damped with increasing |z—z’|. Such a scalar Green’s function is given 
by © 








Y lee! 1 , , | 4) 
Arle’) = Ys -dnltoybnla’sy' Jeers, (8) 
n=1 ” 
where the ¢,,’s are the normalized eigensolutions of the differential equa- 
tion 
hn , © a 
ex? +3 + ke 2¢3 = a 0 


with zero boundary conditions. a k,,’s are the eigenvalues. 

To obtain the integral equations for E, we take the scalar product of 
equation (2a) with N, = N.a, where a is an arbitrary vector introduced 
in order to simplify the handling of the dyadic Green’s function, and 
subtract from the result the scalar product of equation (3) with E from 
the left and a from the right. Integration of the result over the volume 
between the reference planes, chosen far enough away from the obstacle 
so that all the cut-off modes excited by the obstacle are of negligible 
amplitude, yields 


E(r) = Eyne(t) +@%p9 | N47 (r'|r).€' Er’) de’— 


—jw [ (V'xN(r'\r)]?.[w'-H(r’)] dr’, (9) 
where E,,, is the electric field of the wave incident on the obstacle, and 
N? and (YN)? are the transpose of N and (YxN) respectively. To 
obtain the magnetic field, we may utilize the above equation and Max- 
well’s equations, or repeat the above-mentioned steps utilizing equation 
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(2b) and the magnetic dyadic Green’s function M = M,+M,. The 
non-divergent part is 
a HY (p)HP(6')e-™) ( > 
M,(p, 2\e 2)= z B,,X (1)/ a’ \os ~2’ ’ 
n=1 | H2(e)HD(p’)je™-) (zx <2’), 


k2/k? { H2 dS (for TE modes) 
where B, = — — «| : 
“Yn | k2,/we? [ E2?dS (for TM modes). 
] 
The curl-less part is M, = — VV orlr’), 
- | pa - 
where g(r r ) — 7 2k P(x, Yb, (x Y Je-Kmit—2 ° 
m=1 ™ 


The &’s satisfy the differential equation 


Om ihe bn +k, = 0 
oz? lOy*® en 


and the boundary condition n.Yy,, = 0 at the boundary of the guide, 
since M satisfies the boundary condition n.(V x M) = 0. We find 


H(r) = H,,,(r)+ we, [ M?(r' |r). w’.H(r’) dr’ + 


+ jw f (0'xM(e'|e)}" Le’ E(r’)] de’. (10) 


To obtain explicit expressions for the elements of the scattering matrix 
we take the limit of equation (9) as the observation point approaches the 
reference planes. Near the first reference plane (z - —oo), assuming for 
simplicity that only the first mode propagates within the guide, we obtain, 
utilizing the expanded form of the dyadic Green’s function, 


E(r) + E,,.+4,d, E@e, (11) 
where J, = wpy [ (E\!). ce’. E—H. p’. H)e-" dr’. (12) 
Near the second reference plane (z > +00) we similarly find 

E(r) > E,,.+A, J, EMe-”™, (13) 
where J, = wo [ (E?.e’. E—H®). wp’. He” dr’. (14) 
Setting E = E"(r) 
when Eine = Ee-*™, 
and E(r) = E®(r) 
when Eine = E@eti, 
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we find near the first plane that 
E®(r) = EWe-*+ A, J, Ee 
and E%(r) = E@e*+ A, J, E@etir | 
where Ji, = wg [ (Ei). €’. EX—H. p’.H®)e- dr’ 


. 


and Jy = wg | (Ef?.€’. B® HY. p’ He“ dr’, 


Similarly, near the second plane we obtain 

E%(r) = Ee 4A, Jy ood 
and E®(r) = E{Pe*+- A, Jy, Ee 
where Jy, = wt f (E®. e’. EX—H®, pw’. Hei dr’ 


and Jo. = wt [ (E®.¢’. E®—H®. p’.H®)em dr’ 
or ; 
Jum = Wt ( [EM .€ .EM™—H®. p’ HJ" dr’ (n,m = 1,2). 
(20) 
The effect of the obstacle on the propagating mode of the guide is thus 


given in terms of the four J,,,, which are related to the elements of the 
scattering matrix. We note that 


Sy = Ayn =" 

Si. = 14+A,d,. = 1+, 

Sy, = 1+A,Jy = 144, 
and Sy, = Ady = f 


, (21) 


where the r’s and t’s are the reflection and transmission coefficients of 
the obstacle. The subscripts (1,2) on the reflection coefficients r and the 
transmission coefficients t, indicate whether the incident wave is coming 
from the left, z = —oo, or the right, z = 00. The exact solution of integral 
equations (9) and (10) satisfied by the electromagnetic field required for 
the evaluation of the elements of the scattering matrix, is a formidable 
task. Fortunately, however, it is possible to construct stationary expres- 
sions for the elements of the J matrix, thus enabling us to obtain good 
approximate results for these elements whenever a good trial function 
suggests itself. 


3. Variational principle for J,,, 

To construct a variational principle from which integral equations (9) 
and (10) for the electromagnetic field within the obstacle are derivable, 
it is necessary to define an adjoint electromagnetic field. This field we 
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choose as the solution of the differential equation 
VXE! = joH'.p | 
and VxHt = —jwEt.e) 
and assume the adjoint field to have an e~/ time dependence. 
In some problems we shall find the adjoint fields, E' and H', to be 


simply related to the fields E and H. For example, in the case of non- 
dissipative obstacles, where € and p are hermitian (3), 


Et= E* and H' = H*. 


(22) 


For lossy ferrites, if we take Lax’s (4) form of the tensor magnetic 
permeability, having its diagonal elements even in w, and its off-diagonal 
elements odd in w, 


H'(r,w) = H(r,—w) and Et(r,w) = E(r, —w). 


To solve for the adjoint electric field, we introduce the adjoint electric 
dyadic Green’s function, N‘(r’|r), which we find is equal to 


N(rir’)? = N(r’ |r). 
Analogous to finding the integral equation for E(r), we obtain 
Et(r’) = Ejyc(t’) +079 | Et.e’.N(r’|r)? dr + 
+jw [ H.p’.[V xN(r\r’)] dr (23) 
and 
H'(r’) = H},.(r’) + we, | (Ht. w’).M(r’ |r)? dr — 
—jw { (Et.e’).[Vx M(rir’)] dr, (24) 
having used the reciprocity relation 
VxN(r\r’) = [V'x M(r’|r)]’. 


Once again it is necessary to distinguish the solutions for different 
incoming waves. We should note, however, that for the adjoint case, 


Ee and —He-in- 
are the fields for waves incident from the right (z = 00), while 
Ee and —HPeir- 


are the fields for waves incident from the left ( = —oo). Near the first 
reference plane we thus find 


+ , a. (1) o—jyi2’ | + —jyi2' 
E® (r’) = E} e-iy1 +A). 1. Ee iy 


< (2)t i. 2) piy12' 1 + —jyie’ 
and E®(r’) = Ee +A, J}, BMe-m 
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while near the second reference plane 

E(r’) = Ee" +A, J], E@ein” 
and Et (r’) = EPeiv' +A, J}, Ee, 
where 

Thin = wi%ttg [E™t.€ EM +H, p! HM Jee dr. 
Utilizing integral equations (9), (10), (23), and (24) we can show that 
J, = Jt. 


To show this, multiply equation (9) by (E™*.e’) and equation (10) by 
(H™*. w’), to obtain 


[ Ete! EM dr = [ Em. e' Epel-tine dr + 
+w ff Em .¢’.N7(r’ |r). €’. E™ drdr’ — 


—ju | J Ete’ [Vx M(r\r’)]. p’.H™ drde’ 


and 
[Hot pH de = [HOt p! HMw dr + 
+wre, | | Ht, wo’. M2 (r’ |r). p’.H drde’ + 
+jo [ H*. p’. [Ux N(rir’)].€'.E drdr’. 
Similarly from equations (23) and (24) we find 
[ Bm te EM dr = [ Bym.e! Emeline de + 
+w%p f Em ¢’ NT (r’|r). €’. E™ drdr’+ 


+jo fH p(x N(r[r’)].€'.E% drdr’ 


and 
[ Ht. p!.HO dr = — [ Hgm. pw! Hiern dr + 
wie, [ Ht. ’.MM\(r'|r). p’.H ded’ — 


—jo [ Em. e. [Vx M(rir’)]. pH drdr’, 





whence it follows that Jan = Ihn = w'pDan (25) 
Tran J hen 

=Jt — —m ’ 26 

or Jinn Tran an Dan (26) 


5002.44 Ff 
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where 


D,,, [ E™t : e’ é E’™® dr + | Ht : py’ H™ drt wd 


—w pg | f Ete’. N7(r'|r). €'. E™ drdr’— 


nti { (Ht. py’. M7 (e' |r). p!H™ dede’ + 


asp 


+jw [ [Ee [0x M(r[r’)]. pH ded’ — 
—jw [{ H+, pw’ .[UxN(rir’)].€.E™ drdr’. (27) 


It is readily shown that the first-order variation of the right-hand side 
of equation (26) with respect to the fields, subject to condition (25), yields 
the integral equations for the electromagnetic fields within the obstacle. 
Since expression (26) is amplitude independent, condition (25) is auto- 
matically satisfied. Equation (26) is therefore a variational expression 


for the approximate evaluation of the elements J,,,,,. 


4. Application of the variational principle 

The actual application of the variational principle to practical examples 
will be presented in a forthcoming paper. In this section we shall indicate 
some of the different ways in which equation (26) can be utilized in order 
to obtain practical approximate results for the reflection and transmission 
coefficients of an obstacle, 

The full use of the variational method consists of the substitution of 
trial functions with unknown variational parameters into the variational 
expression (26), and the subsequent determination of the values of these 
parameters which make the variational expression an extremum. This 
may seem like a tedious task, and the usefulness of the variational method 
thus be questioned. However, unless one is interested in obtaining only first- 
order results for very small obstacles, this may be the only feasible, in 
any case the most fruitful approach. 

If one is interested in obtaining only first-order results, one may do so by 
simply substituting completely determined trial functicns into equation 
(26). Alternatively, equation (20) may be used, choosing a trial field which 
satisfies equations (25). Since D,,,, is proportional to the square and J,,,, 
to the first power of the amplitude of the field, we can always find an 
amplitude which will make the trial field satisfy equattons (25). For a 
given trial field, the one whose amplitude is so adjusted will yield the 
most accurate results. 

Equations (20) are exact formulae, which have been used by Nikol’skii 
(5) to compute the phase shifts produced by small (obstacle dimension 
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< wavelength A) gyrotropic obstacles in waveguides utilizing trial func- 
tions which he obtains from a quasi-stationary approximation. This 
approximation restricts his results to lossless media (5) and makes his 
result dependent upon the amplitude of the trial field. 

A seemingly different approach has been employed by Berk and Epstein 
(6) to obtain first-order results for the problem of a ferrite post in a 
rectangular waveguide (Fig. 1). Their result, however, is another example 


———» 





X= 77 


FERRITE POST 


/ 





Z=0 
Fic. 1 


of the application of equations (20) with a slight modification. This 
modification involves changing the volume integrals in equations (20) to 
surface integrals. 
The wave equation satisfied by the propagating electric field 
Ei) = E()e-V'ine 
is VxVxEP—FrEP = 0. 
Multiplying this equation by E, equation (2a) by E{?, and subtracting, 
we find 
V .[E® x (Vx EP) — EP x (V x E%)] 
= wi) EY. € . EY —jwEP.V x (p’.H) 
= wp EP. €’. EY—jw(V x El). (pe .H®)+ 
+jwV .[EP x (p’.H)]. 
Integrating over the volume of the obstacle and setting 
VxEP = —jop)Hy 
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yields with the help of the divergence theorem 
wo] [ EW. €’. EY dr — H(). (u’.H®) ar] 
= —Jopo | [E® x HY —E? x H®)].n dS. (28) 


One may thus evaluate the reflection and transmission coefficients 
from the tangential components of E”) and H® on the boundary of the 
obstacle by using equations (20) as perturbation formulae. Berk and 
Epstein used for their approximate field, the long wavelength limit of 
the tangential components of the free space scattering problem. 

To show that the surface integrals appearing in equation (28) are indeed 
the same as the ones used by Berk and Epstein, we note that Berk and 
Epstein restricted themselves to the problem of the scattering of the 
lowest TE mode, having an electric vector parallel to the rod axis 
(y-axis) and independent of the y-coordinate. The symmetry of the 
problem indicates that E% likewise points in the y-direction and is 
independent of the y-coordinate. Furthermore, since the tangential com- 
ponents of EY) and H” are continuous, we may utilize the external fields 


for which , 
—jwpoHQ, = Vx EQ. 


For y-independent fields pointing in the y-direction 


, ' eE®, 
— japon x HQ) = nx Vx BQ, = —— 
on 
on the curved boundary. Thus 
» OOP 4) ED) 1 P 


where & is the magnitude of E, and the integration is performed only 
over the curved boundary of the rod as n x E is zero at the guide boundary. 
If like Berk and Epstein we now set 


6) = PY +6), 


afi (3) 
we obtain J = [(ev Seep oe) dsdy (30) 
J or 7 ae 
» ob of@ _.. 
as QTE _ CF gp) dedy = 0. 
{( * én on ef) - 


The y-integration may be omitted if it is likewise omitted in the evaluation 
of A,. Formula (30) is the result one obtains from formula (6.6) of Berk 
and Epstein if only the contribution from the first term is used in the 
eigenfunction expansion of the Green’s function, as was done by Berk 
and Epstein in their section 8. Again we note that the results of Berk 
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and Epstein are restricted to very thin rods not too close to the guide wall 
and are dependent upon the amplitude of the trial field. 
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SUMMARY 

A helical coordinate system is described and its elementary properties discussed. 
Several electromagnetic problems are solved as illustrations of the use of the system. 
Its advantages are twofold: firstly, it enables problems involving helical symmetry 
to be solved exactly, and secondly, it shows in a new light problems which have 
previously been solved in cylindrical coordinates. The closing section discusses 
points arising out of the examples treated ; in particular it is found that the charac- 
teristic equation for wave propagation in a helical system is always identical with 
that for an appropriate cylindrical analogue, if the usual condition that fields must 
repeat at intervals of 27 in the azimuthal dimension is removed. 


1. List of principal symbols 

Tue following list defines, with their usual meanings, the symbols most 
commonly used throughout the paper. Other symbols which occur only 
occasionally, or any of the following which are used with meanings different 
from those given below, are defined in the text where they occur. Where 
the meaning of a symbol is obvious it is not defined. 


= 


, 6,2 eylindrical coordinates. 
,¢ helical coordinates. 


s 


pitch of the helical coordinate system. 

pitch angle of the helical coordinate system, defined by 
tana = p/2zr,0 <a< Ler. 

parameter giving the z or { dependence of electromagnetic 
fields. 


re 'S 


cee) 


n parameter giving the ¢ dependence of fields. 

q parameter giving the 0 dependence of fields. 

v an integer giving the { dependence in certain cases. 

a radius of a cylindrical surface; when two are present, the outer 
one. 

b radius of helical wire or an inner cylindrical surface. 


c ¢-dimension of helical waveguide, and length of air space in 
helically supported coaxial line. 
{Quart. Journ. Mech. and Applied Math., Vol. XI, Pt. 4, 1958] 
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electric field. 


magnetic field. 
arbitrary vector. 
arbitrary scalar. 
permittivity and permeability, respectively, of free space. 
free-space wavelength, 
= weg ug Br, 
X,(x) linear function of the Bessel functions J,(x) and Y,(z). 


dX, 
‘dx’ 


Xi(2) = 


2. Introduction 


The helical coordinate system treated in this paper, and which the 
author believes has not been previously discussed, is simply related to 
the cylindrical polar system, and may in fact be applied to any problem 
to which cylindrical polar coordinates are appropriate. An example of 
this is the cylindrical waveguide, which is normally treated in cylindrical 
coordinates. The so-called normal modes are in fact degenerate, each 
being the result of a superposition of two helical modes. This is, of course, 
well known, but the helical modes are usually spoken of as ‘circularly 
polarized’. It will be clear from the treatment given below (section 4.2) 
that a more proper term would be ‘helically polarized’. 

An example of a problem to which the present helical coordinate system 
is not suited is the twisted waveguide, which Lewin has treated in a helical 
system related to the Cartesian system (1). It might therefore help to 
distinguish the two systems if we call Lewin’s coordinates the ‘rectangular 
helical system’ and those discussed in the present paper the ‘polar helical 
system’. We shall, however, not refer again to Lewin’s system in this 
paper, so there will be no confusion if we drop the term ‘polar’. 

We shall first give the fundamental theory of the helical coordinate 
system, then solve several electromagnetic problems in the system, and 
suggest other problems, not only electromagnetic, to which it may be 
appropriate. Finally, we shall discuss various points arising out of the 
examples treated here and prove a theorem which enables a great simpli- 
fication to be made in a helical problem by reducing it to an equivalent 
cylindrical one. 


3. The helical coordinate system 

3.1. Description of the system 

Consider first a set of cylindrical coordinates r, 6, z (Fig. 1). The coordi- 
nate surfaces, which are orthogonal, are a set of planes intersecting in a 
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common axis, the z-axis, a set of parallel planes, orthogonal to the first 
set, and a set of concentric cylinders with the z-axis as axis. These families 
of surfaces are given respectively by @ = constant, z = constant, and 
r = constant. 

In the helical system we shall retain the family of cylinders r = constant, 
with meaning unchanged. In place of the parallel planes z = constant, we 
shall take a family of helical surfaces given by z = constant+p0/27, p 


nN 























c | 
B 
Pd 
, a 
a) 
Yl 
5a _—Reference surface for 
the helical system. 
P ee sod 
te 
~ 
Pe a 
y Z Ye 
4 4 
w” Reference surface for the 
cylindrical system. 
om 
0 
Fic. 1. Helical reference surface for right-handed coordinates. 


being the pitch; the surface z = p@/2z7 is illustrated in Fig. 1. The coordi- 
nate z in the cylindrical system is replaced by [ = z—p0/2z, i.e. { is 
measured parallel to z, from the surface z = p@/2z, which in the helical 
system becomes £ = 0. Our third set of coordinate surfaces is the set of 
planes @ = constant; however, we shall no longer use 6, but ¢. ¢ is 
numerically equal to 6, but whereas @ is measured in a plane z = constant, 
¢ is measured in a helical surface £ = constant. It will be seen that this 
helical system is not orthogonal. 

One further quantity that will be useful is the pitch angle a, which is 
the angle made with a plane z = constant by the tangent to the curve 
¢ = constant, r = constant. Its value is given by equation (3) below. 

It might be thought that it would be more convenient to use the system 
consisting of the surfaces r = constant, and ¢ = constant, described above, 
together with a third set of helical surfaces 4 = constant, the % surfaces 
being orthogonal to the ¢ and r surfaces, and left-handed if the ¢ surfaces 
are right-handed. However, the analysis in such a system is very difficult; 
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in particular the wave equation is complicated by the fact that the reso- 
lutes of V?V are not the same as those of the operator (div grad)V. 


3.2. Helical spaces 

Geometrically, we may reach a point (r,6,z+-p) from the point (r, 0,z) 
in one of two ways; we may travel along the straight line r = r, ¢ = 8, 
i.e. suffer a translation in the ¢-direction; or we may travel along the 
helical path r = r, = constant, i.e. suffer a translation in ¢. These two 
geometrically equivalent changes of position may or may not be physically 
equivalent, according to what materials may be present, in what positions, 
in the system under consideration. We may think of the equivalence or 
non-equivalence of the translations as dependent not on the materials 
present but on the properties of the space we are dealing with. This enables 
us to express boundary conditions in a convenient way; according to the 
problem under consideration, we choose the appropriate type of space, 
and the geometry then takes care of itself. 

We shall now define our types of space. If the two translations described 
above are equivalent, we shall describe the space as multiply connected. 
If the translations are not equivalent, we shall describe the space as simply 
connected. There is a third type of space; its definition will be better 
understood after some discussion of simply and multiply connected spaces. 

In simply connected space, the helical surface { = 0 constitutes a barrier 
which cannot be crossed; it is impossible to reach B from A (Fig. 1) by 
travelling parallel to the z-axis, but it can be reached by varying ¢. 
A point can be expressed in one way only, by suitable values of r and ¢, 
and a value of { lying between 0 and p. 

In multiply connected space, the helical surface { = 0 is not a barrier, 
and B can be reached from A by changing { an amount p as well as by 
changing ¢ by 27. If A is the point (r,¢,¢), B may be represented by 
(r,6+ 27,0) or by (r,¢,¢+p). Thus any point can be represented in an 
infinite number of ways by a unique value of r and a value of { which 
depends on the value assigned to ¢. 

The regions between the surfaces [ = np and { = (n+-1)p, for all n, are 
contiguous in multiply connected space; in simply connected space they 
are not, so that by crossing the surfaces { = 0 or { = p, we go outside 
the system. In multiply connected space, when we cross the surfaces 
¢ = Oor { = p, we move to another part of the system. 

The third type of space is that in which { is not limited to the range 
0 <{<>p, but a change of p in { is not equivalent to a change of 27 
in ¢. This type of space may be called disconnected space. It is probably 
of no physical interest, and will not be considered further in this paper. 
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We shall treat below several electromagnetic problems; the coordinate 
surfaces will be so chosen that surfaces of dielectric and conducting media 
are expressible as a single coordinate equated to a constant, or so as to 
coincide with equipotential surfaces. We shall postpone further discus- 
sion of types of space till section 6, since the arguments will be easier to 
follow in the light of the examples. 


Left- and right-handedness 

When referring in future to cylindrical coordinates, these will be taken 
to be right-handed—that is, on looking along the z-axis in the positive 
direction, @ is measured in the clockwise sense. We shall measure { and 
¢ in the same sense as z and @. If the helical surfaces are right-handed, 
as in Fig. 1, we shall, in making a positive change in ¢ with constant (, 
make a positive change in z, so that the pitch p is positive. Similarly, if 
the helical coordinate surfaces are left-handed, p is negative. 


3.3. Elementary geometrical properties 
The following formulae may be derived fairly easily by elementary 
methods, and we state them here without proof. 


Coordinates 
Although these have been fully described in section 3.1, we give the 
equations again here for reference. They are, for the helical system in 
terms of the cylindrical system, 
r=f 


¢ = @ in magnitude, 


but is measured along a different curve, namely, }- (1) 
6 
{a2 
“7T 


For the cylindrical system in terms of the helical system, we have 


2= 1420 (2) 


1 on 


We also have the pitch angle a, given by 


Pp 
tana =5— (0 <a< }p). (3) 

Line elements ds, = dr 
ds, = 4/(r?+-p?/4n*) dd }. (4) 


ds, = dt 
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Surface elements 

dS, = m.ny(r?-+-p*/4n*) ddd = rd 

dS, = drdf 

dS, = \(r?+-p?/4n®) dddr 
Here m and n are unit vectors parallel to ds, and ds, respectively. 
Volume element dV = rdrdddt. 


Distance element 
de? = dr®+-(r2-4-p*/4®) dd? +-dl2+ 2(p/2n) ddl. 

Differential relations 7 
a _ oe 


=—=1 
dz 


> 


$Y | 
Se Ale 


@ 


| 
>! 


le S| 


_ 0 a 8 a 


oe a. we 


51 eG] 
oo 


a 
=| 8 





a 
26 


3.4. Vector relations 
Vector components 

Since the radial component of a vector is perpendicular to both @ and z 
and to ¢ and {, it will be the same in helical coordinates as in cylindrical 
coordinates. Let us consider now a vector V having components V,, V,, 
and confine our attention to a surface r = constant. The components of 
V in both the helical and cylindrical systems are illustrated in Fig. 2. It 
may readily be seen that 


(9) 


Vz = Voseca; Vg = Vg cosa 
V,=V,—Votana; V,=Vi+Vysina 
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Gradient 
In a cylindrical coordinate system, we have 


(grad), = a 














or’ 
1 ap 
hdl, we oor. 
(grad y)g > 20 
ous 
(grad ys), = a 
v.4 Vz 
V 
é 
& > 
Vo 


Fic. 2. Vector components in a surface of constant r. 


If we let grady be the vector V, we can calculate its value in helical 
coordinates by the use of the relations (9). We thus have 


(grady), = © 


(grad b)4 = ah 


rcosa 00” 
tan « dy 
inn ee 
(grad 4) a. —— ie 


Substituting for @ys/20 and @/éz from equations (8), we obtain 


gedit att. } op “) 4-0/4 ‘ene 


or ' cosa ae at 7 4 r oad 


where I, m, and n are unit vectors parallel to ds,, ds., ds. 


+t ana (10) 





Divergence 
We can write 


re 
VV =~ = (h)+(Vgt+VV+V) 


V, and V; can be obtained from equation (10); it is important to remember 
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in forming the scalar product that V,.V; and V;.V, are not zero, since 
the ¢ and ¢ coordinates are not orthogonal. We thus obtain 


osa Vz av; 
cos a $40 : 


divV =~ © < (nV) +2 ag 


(11) 
the cross terms having cancelled out. 


Curl 


The simplest way to treat the curl of a vector V is to write down the 
expression for curl V in cylindrical coordinates and then transform to 
helical coordinates by means of the relations (8), simultaneously replacing 
the cylindrical components of V by helical components according to 
equations (9). The results are 


curl, V = (- oa fe Ww, 1 a) 


rib nr a r op cosa 
(oh 


ato ay —F(%sina)| (12) 


or 





é sin*a &V, lv, p wy 
\./(r?-+-p?/ 4?) ap 6) + ea or -| 


rad 2nr Or 


3.5. The Laplacian operator 
If in equation (11) we replace the vector V by gradi as given by 
equation (10), we obtain 
vy = poh 1 Oy 1 @ 2sina My (13) 
ror "or r? ad? cos*a 80? r cosa alad 


which may be written 


eys\ 1 a p? \e%y ary 

2 je ee A, a 14 
i r wal’: +73 ait ( + Dials (2 alane: (14) 
The symmetry of a physical system will often impose the condition 
2/éb = 0, i.e. the helical coordinate surfaces may be chosen to be equi- 
potential surfaces. In this case, the Laplacian reduces to the simpler 


form ss 2 Sf )+(1 ~ Galo (15) 


The wave equation 


We have (v=, va ~ (16) 
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where v is the velocity of the waves. Let us assume a time variation of 
the form e/, and write w?/v? = ,?. The wave equation then becomes 


dr? por | ridge? \ ' 4n2rtsat? \2nrtlated 





and we shall look for solutions of the form 
fb = R(r) P(d)Z(Q). (18) 

Substituting for 4 from equation (18) into equation (17), we obtain 

, -_ Pp" p? Zz” pP’Z' 

—+—+-—_ 4 + ———. j-—— — -_ 2 — (), 19 

R'?R ) ap+| : caly ar PZt * (19) 
where the primes denote differentiation with respect to r, ¢, or ( as the 
case may be. We cannot now separate the variables, but we can obtain 
an equation in r alone if P’/P, '’Z'/PZ, and Z"/Z are all constant. By 
analogy with the case of cylindrical coordinates, we assume 








4 = et iB 
which gives Zz" 
5 =P) 
* a ae (21) 
PZ’ 
ema 
and equation (19) then reduces to 
| eet . 
a tapt®- a = 0, (22) 
where q=n bs (23) 
and = x*—f?. (24) 


Equation (22) becomes Bessel’s equation of order q if we take kr as the 
independent variable. The solution is thus 

R = X,(kr), (25) 
where X,(x) denotes a linear combination of J,(x) and Y,(z). 

There is always the possibility, of course, that other solutions of the 
wave equation exist which are not obtained in this way, but this applies 
equally to all coordinate systems, since the procedure we have followed 
here is the one commonly used in other coordinate systems. 

The solution of the wave equation in helical coordinates is thus seen to 
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be related to its solution in cylindrical coordinates. The differences are 
that the order of the Bessel function is not the same as the n in the ¢ 
dependence, and n and q are not constrained, in general, to be integers. 
It will often be convenient to choose the helical coordinate surfaces ({ = 
constant) to be equipotential surfaces; in this case, @/é¢ = 0, Le. n = 0 
and q = —p§/2rn. 

It should be noted that in a right-handed system we must take the — 
signs in both equations (20) for a forward wave, or the + signs in both 
equations for a backward wave. We have also pointed out above that p 
is positive. In a left-handed system, we must take one + sign and one 
— sign in equations (20), and p is negative. Thus p8 and n have the 
same sign in either system for forward waves and for backward waves. 
In changing from a wave in one direction to a wave in the other, or from 
a systen. of one hand to a system of the other hand, q changes sign, but 
retains the same numerical value, and R is unchanged, as we should 
expect. The situation may be summed up thus: 














Forward wave Backward wave 
Right-handed Left-handed Right-handed Left-handed 
coordinates coordinates coordinates coordinates 
Pp tive —ive +ive —ive 
B ive —ive +ive + ive 
n ive +ive +ive —ive 
q +ive or —ive —ive or +ive —ive or +ive +ive or —ive 




















The sign of q will depend on the relative magnitudes of n and p§/27. If 
n is less than p§/2z, q is to have the sign given first in each coitumn of 
the table. This will usually be the case, and is so in particular when n = 0, 


4. Electromagnetic theory 

4.1. Fundamentals 

In section 3.5 we have solved the wave equation for an undefined scalar 
independent variable %. In electromagnetic theory this variable will be 
the magnitude of an electric or magnetic field component, and in cylin- 
drical coordinates it would be EZ, or H,. In helical coordinates the variable 
is E;+H,sina or H;+H,zsina. One of these can be put equal to 0 in 
most cases to give H-waves or E-waves. We shall also have x? = w€, fio. 
By virtue of equations (20), (23), (24), and (25), we have, for forward 
E-waves in a right-handed system, 


E,+E,sina = X reins eae 


: (26) 
H;+Hg sina = 0 
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while for H-waves 
E+ Ey sin « = 0 \. 27) 
H,+Hgsina = X,(kr)ePe-indeiot | 
It is possible to express all the field components in terms of EF, and H, 
by separating them out in turn from Maxwell’s equations. This, however, 
is a tedious task in helical coordinates, and it is simpler to take the ex- 
pressions for the field components in cylindrical coordinates and then 
convert to helical coordinates with the aid of equations (8) and (9). For 
a system in vacuo we obtain 


1 fap q : .p[OHr . 2 " 
E. = Bl — [H+ Hysin a] +j6| 8+ 5 (Bgsin || 


a ; Bq 
—(H,s —=[E,+ Es 
anf gsina)| na a sin al] 
OH; 


C WE 7 —— 
s+ = - (Hy sin ~)|- ~~ [By + Eg sin 5} 
H, = ee |. Pa H+ Hi, sin *) jes! on 5; (Hg sin «|| 


k?cosal r 





Substituting from equations (26) and (27) into equation (28), we obtain 
for forward E-waves in a right-handed system 


EB. = * (i (ker)e-1BSe—ine 


A X ,(kr)e*Pe-ind 


kr 
By = (1+ 74 ) X ,(kr)e~tBle-ind 


H, = 504 X (kr)e-iPle-im4 


H, = = poets Xi (ker)e-HBle-ing 





= jwe Py le —j —jn 
B= SE rie et 
where X'(kr) = dX,/d(kr), and the factor e/ is understood; for H-waves 
we have 


Z, = “ee zx qikrye e-iBle—ind 


_ Jp ~iplp—ing 
= YY qikr)e- e~ 





_ —IJMOP ¥ (r\o-1Ble-ing 
ob X fhr)e—Pte 
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H, =P x (er)e-iPhe-ind 


H, = ——P4_ x (br\e-iBle-ind 
$ rence” , 





I on ( + Peal X lerje-iBhe-ind | 
For a left-handed system, replace p by —p and n by —n in equations (29) 
and (30). For backward waves, replace 8 by —f and n by —n. For 
backward waves in a left-handed system, make both these changes; the 
net effect is to replace p and B by —p and —f respectively, the two 
changes in n having cancelled. In making these changes, it must be re- 
membered that q is a function of p, 8, and n. Equations (29) and (30) apply 
equally to a system filled with a homogeneous isotropic dielectric or mag- 
netic material, if «, and py are replaced by ee, or po. 

We shall now treat various examples; in the main, the problems amount 
to substituting the appropriate boundary conditions into equations (29) 
and (30) and so deriving a characteristic equation. However, it is not 
always 8 that the characteristic equation is to be solved for. 


4.2. The cylindrical waveguide 

Since the system contains the axis r = 0, at which the fields must remain 
finite, the solution of the wave equation that involves Y,(kr) is excluded, 
so that the field equations are as given in equations (29) and (30) with 
X (kr) replaced by J,(kr). The helical coordinate surfaces are conveniently 
chosen to be equipotential surfaces, so that n = 0 and p is thus equal to 
q wavelengths. Now 

e-iBl — e-iBtg—ipBbi2n — ¢-ipzg-iad 


so g must be an integer and p is therefore an integral number of wave- 
lengths. 

We now apply the boundary condition that EZ, = 0 at r = a, a being 
the radius of the waveguide. From equations (29) and (30) we thus obtain 


Jka) = 0 (31) 


for E-waves, and Jka) = 0 (32) 


for H-waves. These are the characteristic equations, to be solved for k 

and hence 8. They are identical with the characteristic equations obtained 

by the treatment in cylindrical coordinates, but do not refer to the same 

kinds of waves. In cylindrical coordinates, we should have been dealing 

with waves in which a given point in the wavefront moves parallel to the 

axis of the guide; in the present case, a given point in the wave front 
5002 .44 Gg 
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moves along a helical path. Such waves are commonly said to be circularly 
polarized; it is clear, however, that a better term would be ‘helically 
polarized’. Circular polarization only exists in an infinite homogeneous 
medium. 

In a cylindrical waveguide, there is no reason why one sense of helical 
polarization should be preferred to the other. Let us therefore consider 
the case of two oppositely helically polarized waves, of the same frequency, 
having equal amplitudes, travelling together in the same direction along 
a guide. We denote by S any one of the field components, for an Z- or 
H-wave, and write for the right-handed wave 

S, = Kf(rje®, 
which may be written 
S, = Kf(r)e~iBe-i09. 


For the left-handed wave, 


S_ = 4K f (reel, 
where the + or — sign is to be chosen according to the values of the 


constant K and of the radial function f and the parity of g. For the 
composite wave S,+-S_, we obtain 


S = Kf(r)e-i8=(ei49 + ¢-i0), 


ie. S = Kf(r)e-iBel°*\a0, (33) 
| sin} 

which is identical with the field component found by the treatment in 
cylindrical coordinates. This demonstrates the degeneracy of the normal 
modes of a cylindrical waveguide, each of which is a combination of two 
helically-polarized waves. Apart from this aspect, the example is, perhaps, 
trivial; it does serve, however, to give us some feeling for the helical 
coordinate system before going on to some more difficult and interesting 
problems. For this case, it can be seen that the helical-coordinate treat- 
ment is no more difficult than the cylindrical-coordinate treatment, and 
this suggests that the former may be applied without loss, and perhaps 
with some gain, to problems previously treated in cylindrical coordinates. 

4.3. The helical wire 

We shall now consider the problem of wave propagation along a helical 
wire in otherwise free space. The wire is supposed to be infinitesimally 
thin, and we choose the helical coordinate system to have the same pitch 
as the wire. The equations of the wire can then be written 


kee (34) 
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The boundary condition to be satisfied is that the tangential component 
of electric field on the boundary is zero, i.e. 


Ey+E;sina = 0. (35) 
Consider first H-waves. From equations (29), equation (35) becomes 
(1—B?/k*)J,(kb)e-1P%e-in? sin a = 0. (36) 


This condition can be satisfied in one of two ways, either by putting 
J,(kb) = 0 or by adding a forward and a backward wave so as to replace 
eB. with sin(27m{/p), m being an integer. This gives two sets of waves, 
and in general the field configuration will be given by adding members of 
both sets. If the conductor has extension in the surface r = 6, the sine 
solutions vanish and we must choose J,(kb) = 0. If the conductor has 
extension in the surface { = 0, the condition J,(kb) = 0 does not apply, 
and we must choose the sine solutions. 

If a second conductor is present, we have a transmission line; an 
example of this is the slow-wave tube to be considered in section 4.4. 

The case of H-waves is similar, except that the condition J,(kb) = 0 
that we have for E-waves is replaced by J)(kb) = 0. 


4.4. The slow-wave tube 

The idea of the slow-wave tube was first formulated by Kompfner (2, 3). 
Many devices have since been invented for producing slow waves, but in 
Kompfner’s form the slow-wave tube consists essentially of a helical wire 
of radius 6 and pitch p, concentric with a tube of radius a, a being greater 
than 6. Electromagnetic waves may be thought of as travelling along the 
. wire with phase velocity vg, so that if the pitch angle of the wire is a) the 
component of velocity in the axial direction is 

Vg = Up SiN ay. (37) 

Kompfner found experimentally that vy is approximately the velocity of 
light (3). 

Since we have two distinct conductors, the device may be treated as 
a transmission line; it is possible, of course, to treat it as a waveguide, 
but as a slow-wave tube it functions as a line, and we shall only consider 
this case here. We shall proceed by analogy with the treatment of the 
coaxial line given by Lamont (4), and for simplicity assume that the wire 
is infinitesimally thin. 

Anticipating Lamont, we assume a two-dimensional potential function 
ys, of the form ob = (rend, (38) 
For the electric field, we have 


E = —grady, (39) 











R. A. WALDRON 


" Op ing 


i 


or 


E, = J ¢-ing 


Meno 
rT COS a 





—jntan a 
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The magnetic field components could be obtained from Maxwell’s equa- 
tions by substitution from equation (40), but this is rather tedious and 
we do not require the resu'ts for our present purpose. 
We require div E = 0. Substituting from equation (40) into equation 
(11), with V replaced by E, and equating to zero, we obtain 
ay Leap n? 
6 rer 


b= 0. (41) 


The solution of this is ys = Ar*”-+-constant, 
and since % must remain finite at r = 0, we reject the negative sign. Also, 
the boundary conditions must be satisfied that Ey and E; vanish at r = a, 
which enables us to put 4 = 0 at r = a, ie. the potential of the outer 
conductor is zero. Thus we write the solution of equation (41) in the form 
b = a"—r", (42) 
The function e~/#* must also satisfy the wave equation for propagation 
in the z-direction. Now, 


ow whe jBe (rye ings iBle—ipBdl2m 


Oy Lab (n+pB/27)? 3 2 B wo 
Sp MBO pf lo lar = 0 
(43) 
where v, is the velocity of electromagnetic waves in the z-direction. Sub- 
tracting equation (41) from equation (43), we obtain the characteristic 
equation £2 = w?/v?, 
i.e. v, = w/B (44) 
as we should expect. 
Here f is fixed by the requirement that at ¢ = 0 and { = p the fields 
take the same values, so that p8/27 must be an integer, say v, and we have 
v, = wp/2nv, (45) 
so that v, is fixed by the frequency and the pitch of the helix. Equation 
(45) can only be expected to hold as long as p is less than vA). 
From equation (45), bearing in mind that w/(€yp5) = w/e = 27/Ap, 
where c is the velocity of light, we obtain 


v, = pc/vaAg, (46) 
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which shows that slow waves exist when p < vAy. The velocity v, may be 
regarded as the z component of the velocity of waves travelling along the 
wire. We have p 
siIN a = —— >> 
X a (427°b? +- pp?) 
(4722+ 2 
so that ee (47) 
vg 
Thus the waves travel along the wire with the velocity of light if the length 
of one turn of the helix is equal to vAy. Unfortunately, Kompfner’s papers 
do not give enough information to enable us to check whether this condi- 
tion was satisfied in his experiments. 


4.5. The helical waveguide 
This problem has been previously treated by the present author by an 
approximate method. The waveguide is bounded by the surfaces r = b, 
r=a,(= 0, €=c,c being less than p. We shall consider Z-modes and 
H-modes separately. 
E-modes 
In equations (29) we put 
X (kr) = J,(kr)+ AY (kr). (48) 
On €=0 and (=c, £;,+E;sina = 0. Thus we must take a forward 
wave and a backward wave and combine them so that in equations (29) 
e~'B¢ is replaced by sin vzf/c or cosv7f/c as the case may be. We also have 
K? = wey py— B® = weg py—v*x*/c*. (49) 
The surfaces £ = 0 and { = c are equipotential surfaces, so that 2/é¢ = 0, 
i.e. n = 0. The field equations thus become 


vm Jer) +-AY ,(er)jsin 4 


; = et 1) + AY (ker) sin 76 
E, a {J,(kr)+ AY, (kr) }sin : 


E, = (1+ +5728 Bis) dr) + AY,(kr)}oos == 


— WE] 


—— (I(r) + AY, (kr) }e0s 


ate © (Telkr) + AY‘ (kr)}cos . 





= font * (Jelkr) + AY ,(kr)}oos 
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We must also satisfy the condition that (Z,+ #;sin «) vanishes on r = a 
and r = }, i.e. J,(ka)+-AY, (ka) = 0, 
J (kb)+-AY,(kb) = 0. 
These two equations must both give the same value for A, so that 
Jka) Y (ka) 
J(kb) Y,(kb) 


= 0, (51) 








which is the characteristic equation, to be solved for q, since k is fixed by 
equation (49). 
H-modes 


These may be treated in the same way as E-modes. The characteristic 
equation is found to be 





Jka) Yi(ka)| = 0, (52) 
Ji(kb) Yi(kb) 


and the field equations are 


BE, = CHOI 7 kp) 4 — 
Cc 





k2r ** 
E, = Jobo moe kr) +A’'Y ((kr)}sin vat 
 keos c 
7 jw P I’ (kr v7 
E, = or ah JTkr)+A’'Y,;, (kr) sin(” 5) 
(53) 
H, =. cal Jkr)+ A'yY’ alkr)} eos(" =| 
ke c 
7; juTrq f y\ 1 A’V (hk mt 
Hy = Pt filer) + AY (kr)} ma : ] 
Par m me 
H; ( r spa) (kr) T Y(kr)} sin( Cc 


Now, from equations (2) and (23), remembering that » = 0, we may 


write e-iBE — e-iPegipBdl27 — ¢-iBag—iad. 


Thus qg is the number of waves per turn of the helix. The characteristic 
equations (51) and (52) are identical in form with those for the infinite 
circular guide, which is a guide with surfaces r = b, r= a,z=0,z= Cc, 
with the restriction removed that 6@ is limited to the range 0 to 27. Thus 
the number of waves per turn in the helical guide is the same as in the 
infinite circular guide. This conclusion was reached by the author in the 
previous approximate treatment (5), where the method was to treat the 
infinite circular guide exactly and then apply perturbation theory to find 














A HELICAL COORDINATE SYSTEM 455 


the change in the azimuthal phase constant (i.e. g) on shearing into the 
helical form. The change was found to be exactly zero. 

The solution of the characteristic equations has been discussed by the 
author in the previous treatment, where tables are given of kb as a func- 
tion of q, with the ratio a/b as a parameter. In using these tables, it should 
be remembered that the notation is different from that in the present 
paper, and care should be taken to avoid confusion. 


4.6. Coaxial transmission line with helical dielectric support 

The geometry of this system is similar to that of the helical waveguide 
treated in section 4.5; the cylindrical walls r = a, r = 6 are conducting 
surfaces, and the region 0 < { < ¢ contains only air. The regione < ¢ < p 
(p > c), however, instead of being filled with metal as in the helical guide, 
contains a dielectric, which acts as a support for the inner conductor. 
The system may be treated as a waveguide or as a transmission line; the 
transmission line problem has been treated approximately by Griems- 
mann (6). We shall now give exact treatments of both the waveguide and 
the transmission line problems; the waveguide problem is of interest for 
the comparison between its characteristic equation and that of the helical 
waveguide. 
(a) Waveguide treatment 


The field equations take different forms according as we are dealing 
with the air- or dielectric-filled regions. For the air-filled region, we re- 
place 8 by £, and k by ky, in equations (29) and (30), and have 


ke = wey o— Bo. (54) 
In the dielectric-filled region, we replace 8 by f,, k by k,, and e, by €€,, 


and have ES = cw%ag jig e—Pt. (55) 
We shall not write out the field equations in full at this stage. If we 
change the ¢ coordinate but keep ¢ fixed, the phase change on traversing 
a distance p is c8,+-(p—c)8,. If we make the same change in position by 
keeping ¢ fixed and changing ¢ hy 27, the phase change is 27m. Thus we 
— cPo+(p—c)B, = 2n(n—9), (56) 
where q is an integer. We do not necessarily have the same phase at 
f)-+p as at fy, so that n will probably not be zero. Equation (56) may be 
compared with equation (23). 

At €=0 and {= c, the values of E,+E;sin « and E, must be con- 


tinuous, and this can only be the case if k, = k,. Thus from equations 
(54) and (55) B2 = 2-+cw%p prp(e—1). (57) 
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In future we shall write k to mean either k, or k,. From equations (56) 
and (57) we obtain 
By = 2n{(n—g)e+(p—e),|[(n—g)?+ p(e—1)(2c—p)/A3]} (58) 
and so k is a function of (n—q) which we might write as 
k = k(n—@Q). (59) 


Atr=aandr=b, Ey+£, sina must be zero, so that for H-modes we 





must have im ¥,(ka) | = 0, (60) 
JA(kb) Y¥,(kb) 

while for H-modes 
Ji(ka) Yi(ka)| = 0. (61) 








Sikh) Yi(kb) 

Equations (60) and (61) may be compared with equations (51) and (52). 
However, whereas the latter are to be solved for g, equations (60) and (61) 
are to be solved for k, q being assigned according to the mode. From 
equation (59), (n—gq), and therefore n, can be obtained, and we can then 
calculate 8, and £8, from equations (56) and (58). Unattenuated propaga- 
tion will only occur if 8, and f, are real, which, by virtue of equations (58) 
and (59), will impose a restriction on the solutions of equations (60) and 
(61) which are permissible. However, except for g = 0, 8, and f, have 
real parts even in the attenuation regions, so that they are not true stop 
bands, and some propagation can take place. 


(b) T'ransmission line treatment 
As in the case of the slow-wave tube, we may follow a treatment 


analogous to that of Lamont (4) for the coaxial line with homogeneous 
support. We assume a potential function of the form 


b = or, ¢) (62) 
and write E = —grad pe-iRle-ing, (63) 
Le. E. = _ Of ,-iBty-ing \ 
. ér 
—l Op , ‘ 
E, = - a aoe iBle ind 64 
¢ reosa od ‘ (64) 
E,- = tan a ye e-iBle—ing ) 





r 0d 
ys must satisfy Laplace’s equation, so that 
ay lep 1 ep 


PEs AS Md: Ley 65 
or? rer | od? (65) 
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The wave equation is, in air, 


(V2-+ weq p9)(pe/Po*) = 0, 


ie. (V2+- we, pp) (oe *Pole—ivPodi27) — 0, (66) 
If we write out equation (66) in full, and subtract equation (65), we obtain 
By = weg Ho. (67) 
Similarly, for the dielectric region, 
Cre (68) 


These equations for 8, and f, satisfy equation (57), as we should expect. 
Also, we may expect that 4 is everywhere independent of ¢, as it must be 
on the two conducting surfaces r = a and r = b. Thus EZ, and E; are zero, 
so that there is no difficulty in satisfying the boundary condition, which 


al 


is that Z, is continuous at the air-dielectric interface. The solution of 
equation (65) becomes = Alogr+B (69) 


as in the case of the homogeneously supported line. 

The form of equations (67) and (68) shows that the behaviour of this 
line, as a line, is analogous to that of the coaxial line supported at regular 
intervals by annular disks of dielectric material. It is well known that when 
a line is loaded periodically, and when the wavelength of waves on the line 
is comparable with the spacing of the loads, a resonance effect occurs, so 
that the frequency spectrum is divided up into a series of pass-bands and 
stop-bands. This applies to the line with the annular disks, and therefore, 
by analogy, to the line with helical support. The calculation of the posi- 
tions of the bands is easier for the cylindrical case than fur the helical case, 
and in accordance with the theorem to be proved in section 6, we expect 
the results to be identical. It is not our intention to perform this calcula- 
tion here; our purpose is to illustrate the uses of the helical coordinate 
system, not to give an exhaustive discussion of each topic mentioned. 
The theory of stop-bands and pass-bands is fully discussed by Brillouin (7). 


5. Further applications 


The helical coordinate system is applicable to any problem involving 
helical or cylindrical symmetry, although in the cylindrical case it may 
be that no new results are obtained. We shall indicate here a few problems 
which may usefully be discussed in helical coordinates, without, however, 
attempting such discussion. 

A fluid flowing in a cylindrical tube may follow a helical path; this may 
be due to asymmetry in the pipe, causing rotation of one hand to be 
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preferred to that of the other. Other systems are the flow of air behind 
an air-screw, or of water behind a ship’s propeller. It may be helpful to deal 
with screw dislocations in crystals by means of helical coordinates. 

This coordinate system may also be useful in dealing with the motion 
of an electron in a magnetic field. An electron which initially has com- 
ponents of velocity parallel to and perpendicular to a magnetic field 


follows a helical path; this effect is made use of in the helical focusing of 
charged particles. 


Note on toroidal helical coordinates 


[Imagine a circle [ in a plane I], and a line Z in I] which does not 
intersect the circle. Move the circle parallel to LZ with constant velocity, 
and at the same time rotate II about Z with constant angular velocity. 
Let the velocity of the circle parallel to Z be such that while Il undergoes 
a complete revolution the circle moves a distance greater than its diameter. 
The resulting figure is a helical toroid, related to the familiar circular 
toroid in the same way that the helical wire of section 4.3 is related to 
a circular wire. The coordinates appropriate to this figure are (i) the 
radius r measured from the centre of the circle, in the plane of the circle, 
(ii) the angle @ measured round the circle, (iii) the angle 4 measured along 
the helical locus of a point in the circle, as it performs the motions described 
above. It is evident that the coordinate ¢ in this system is identical with 
the coordinate ¢ in the polar helical system with which this paper is 
chiefly concerned. 

The toroidal helical coordinate system is appropriate to such problems 
as the helical waveguide of circular cross-section, which is probably of 
only academic interest, or the helical wire of finite thickness, which is a 
better model for Kompfner’s slow-wave tube (section 4.4) than the helical 
wire of infinitesimal thickness which we have taken in this paper. 

It is outside the scope of this paper to Give an exhaustive discussion of 
toroidal helical coordinates; but it might be suggested here that, by 
analogy with the theorem proved in section 6, it is probably true that, for 
mathematical purposes, a helical toroidal system can be replaced by a 
circular toroidal system and still lead to the same results. Without going 
into the matter thoroughly, it is impossible to prove or disprove this state- 
ment, but intuitively one feels that it ought to be true. If so, then the 
helical wire of finite thickness can be replaced by a series of equally-spaced 
circular wires along a common axis, with the condition removed that on 
changing @ by 27 one returns to one’s starting-point. This represents a 
great simplification of the problem, but the problem of solving the wave 
equation in (circular) toroidal coordinates still remains. 
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6. Discussion 


In section 3.2 we introduced the idea of simply and multiply connected 
spaces. We can pursue the idea a little further now that we have seen the 
results of certain problems. 

In the cylindrical waveguide (section 4.2), a change of p in { or of 27 
in ¢ brings us to a point at which the fields are identical with those at the 
starting-point. The two changes are equivalent, so this is a case of multiply 
connected space. 

In the helical waveguide, the surfaces {= 0 and { = isutute 
barriers; the changes are not equivalent. A change in { of 2 -8 us out- 
side the waveguide; a change in ¢ of 27 leaves us further along the guide, 
but at a point where, in general, the fields are not identical with those at 
our starting-point. This is an example of simply connected space. 

The slow-wave tube and the coaxial line with helical support are rather 
different cases. The effect of changing { by p or ¢ by 27 is to bring us to 
the same point, but it is not an equivalent point, i.e. the fields are not the 
same as at the starting-point. These are examples of multiply connected 
space, and serve to demonstrate that in multiply connected space, the 
geometry of electromagnetic fields need not conform to the geometry of 
the material system. 

From a consideration of the above cases, the following rules emerge: 

(i) If it is possible to find a line parallel to the z-axis which does not 
enter a conducting medium, the system is multiply connected. This does 
not necessarily mean that all lines parallel to the z-axis must satisfy the 
condition; it is sufficient that it holds for at least some. 

(ii) If it is not possible to find a line parallel to the z-axis which does 
not enter a conducting medium, the space is simply connected. 

(iii) All helical spaces that occur in electromagnetic problems are either 
simply connected or multiply connected. 

(iv) In multiply connected space, the periodicity of electromagnetic 
fields is not necessarily simply related to the pitch of the material system 
unless the system has cylindrical symmetry, in which case the pitch of the 
coordinate system can be chosen to give a simple relationship. 

The practical usefulness of deciding what kind of space we are dealing 
with lies in the fact that for multiply connected space q is constrained to 
be an integer (including zero), and the characteristic equation has to be 
solved, effectively, for the wavelength parallel to ¢, i.e. for 8. In simply 
connected space, on the other hand, f is constrained to be 7 times an 
integer, and the characteristic equation has to be solved, effectively, for 
the angular wavelength along the @ dimension, i.e. for g. A considerable 
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simplification of the treatment results if it is realized at the outset which 
of these cases applies to the problem under consideration. 

We have seen that the helical waveguide has the same characteristic 
equations as the infinite circular guide, and that the treatment of the 
cylindrical waveguide in helical coordinates leads to the same charac- 
teristic equations as the treatment in cylindrical coordinates. Similarly, 
the helical wire in the slow-wave tube may be replaced by a series of 
circular wires, of radius 6, concentric with the tube, and with equal spacing 
p from each other. The same characteristic equation will then result as 
in the case treated above (section 4.4). In Griemsmann’s problem (section 
4.6), the helical dielectric may be replaced by a series of annular disks, 
with the same thickness, p—c, as the actual dielectric, and regular spacing 
p. Again the same characteristic equation will result. In the last three 
cases it must be noted that in the cylindrical analogue, the condition that 
all fields repeat at intervals of 27 in @, which usually applies to cylin- 
drically symmetrical systems, does not apply; the cylindrical analogues 
are mathematical devices that cannot exist physically. 

The remarks of the preceding paragraph suggest a result which we shall 
now prove. 


THeEoremM: If a cylindrical analogue of a helical system be made by re- 
placing the helical surfaces € = constant with planes z = constant-+-sp (s an 
integer), and if in the case of simply connected space the requirement is not 
imposed that the conditions shall repeat at intervals of 27 in the azimuthal 
dimension 0, the same characteristic equation applies to the cylindrical system 
as to the helical system. 


To prove this, we write down first the field equations for the cylindrical 


system. They are: 
H l {e (‘ 7  J@€ eF,\ 
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’ ky? YO or r €z\ 26 ) 
The corresponding expressions for the helical system are given by equa- 
tions (28), and if we bear in mind equations (9), where V may be E or H, 


it becomes apparent that any condition on £,, Ey, +E;sina, H,, or 
H+ H;sin « in the helical system must lead to the corresponding restric- 
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tion on E,, Ey, H,, or Hg, respectively, in the analogous cylindrical system. 
Boundary conditions involving E; or H; always involve them in combina- 
tion with E, or Hy as E,+ Egsin « or H,+ Hg sin a, which are identical 
with 2, and H,. Comparing equations (28) and (70), it is apparent that 
boundary conditions in the helical system must lead to the same equations 
in £, and H, as the corresponding boundary conditions in the cylindrical 
system. The theorem is thus proved for electromagnetic waves; an ana- 
logous proof could be given for any other kind of wave motion. 

The application of this theorem leads to a great simplification of the 
problem of deriving the characteristic equation for a helical system, but 
the field components and normal modes are better discussed in helical 
coordinates. However, the most difficult problem is always to derive and 
solve the characteristic equation, so that anything which simplifies this 
process is valuable. 

To sum up this discussion, the best way to treat a helical system is to 
decide first what kind of space is appropriate, and describe the cylindrical 
analogue. The characteristic equation for the analogue may then be 
obtained, and this is the same as the characteristic equation for the actual 
helical system. Having solved the characteristic equation, the field equa- 
tions for the actual system can be obtained by substitution in equations 
(29) and (30), or, more generally, (28). 
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SUMMARY 


This paper deals with the torsion and flexure of a solid cylinder of simply-connected 
boundary cross-section which can be mapped on the ring-space between two con- 
centric circles, the inner one being the transform of an inner slit or cut in the cross- 
section which, however, is not a physical boundary of the cylinder. The solutions 
for the potential functions, torsional constant, and centre of flexure are determined 
in the general case and particular values of these are found for the elliptic cylinder 
and a cylinder having a circular arc as internal cut. 


1. Introduction 


In his book, Some Basic Problems of the Mathematical Theory of Elasticity, 
Muskhelishvili (1) deals, inter alia, with the fundamental boundary value 
problems of two-dimensional elasticity using various methods, and the 
case of the fundamental problem for a continuous ellipse is dealt with 
by means of conformal mapping on to a ring-space and by using infinite 
series solutions. The Cauchy integral methods described by Muskhelishvili 
and used also by Sokolnikoff (2) for the torsion and flexure problem cannot 
be applied to the particular case of the ellipse when one uses the usual 
conformal transformation of the form 


Wit 
t+ 


Z z(t) = k{ 


(K > 0, m > 0), 

t j 

transforming the space inside the ellipse on to the ring-space between two 
circles, since z’(f) = 0 at points inside the circle corresponding to the 
ellipse. The results for the torsion and flexure problem for the ellipse, as 
examples of the fundamental boundary value problem have not been 
given by Muskhelishvili, and Sokolnikoff determines them by means of 
real function theory methods which have been known for some con- 
siderable time. Complex variable methods using conformal transforma- 
tions have been used by Stevenson (3) and Sherman (4) for problems of 
plane stress involving a thin isotropic elliptic plate, and one of the present 


[Quart. Journ, Mech. and Applied Math., Vol. XI, Pt. 4, 1958] 
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authors (Morris) has also developed complex variable methods for two- 
dimensional potential theory problems. The results of Morris (5) for 
Saint-Venant’s torsion and flexure problem were published in three papers 
subsequently referred to as (I), (II), (III), and the paper (II) dealt 
specifically with cylinders of elliptic type for which the transformations 
used had an internal cut joining branch points of the transformation. 

It was in dealing with aeolotropic plate theory problems later, that 
Morris and Livens (6, 7) had to consider again the question of choosing 
potential functions which were non-cyclic in a regio: -ontaining singular 
points of the transformations used for the particular shaped plates. The 
problem was similar to the torsion and flexure problem as dealt with in 
(11) and at first it was thought that the method of that paper could be 
used. However, it was discovered that this method did not yield the 
correct result (verified by reduction to the isotropic case), and a new 
method had to be evolved. This method for the plane stress problem of 
a thin elliptic aeolotropic plate under given edge conditions is given in 
the paper by Morris (7). Morris realized immediately the possibility that 
the results of (11) might be incorrect, and on examining these it was soon 
found that this was the case. The object of this paper is to indicate that 
error and by application of the new method, correct the subsequent results 
of (LI). 

The new method will be shown to be applicable to any boundary 
section which is such that the double-Riemann surface inside the cross- 
section having a cut along a given curve inside the cross-section is trans- 
formable to a ring-space between two concentric circles, which we define 
in the t-plane by {t 1, |t} = e-**, The cut inside the cross-section 
corresponds to the circle |t} = e-* so that with |t} = 1 on the boundary 
the cross-section inside the boundary in the z-plane corresponds to the 
ring-space e~* < |f/ < lin thet-plane. The particular cases of the elliptic 
cross-section and a cross-section having a circular are as internal cut are 


treated. The notation of (II) is followed unless otherwise stated with 
t ¢ it 


2. The new method 


It was stated in (II), section 4, that the regularity of the potential 
functions could be secured provided the gradients of the potentials do 
not become infinite at the singular points of the transformation. This 
condition was secured and a solution for the complex potential function 
(2, determined as 


Qo, = > {21 eh +O_, ert} + 4(1+4 20)22, 
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where ae * ms 
Q,3 t(1+20)f rt $(9-—9-+)—99 +, 
Q2 rl }(1+-20)f, +g ge - Jr) — 99, 
the coefficients f,, g,, etc., depend on the shape of the cross-section, and 
o is Poisson’s ratio. The relationship between this function Q, and the 
potential function x given by Love (8) is 


x = re[4(2+0)z3—Q,], 


and when the appropriate values of f,, g,, ete., for the elliptic cross-section 
are substituted in the above value of Q,, this is found not to agree with 
the value of y given by Love for the ellipse. This indicates the error in 
(II) and the consequent incorrectness of the results in (I1). 

The required condition is that the value of the potential must remain 
unchanged in crossing the cut joining the branch points of the transforma- 
tion. From the papers by Livens and Morris (6), and Morris (7) it is 
noted that if a characteristic boundary condition for the potential func- 


tion Q(t) is of the form > B,,t%} (where t = t, = e-€ on the boundary), 
* 
0 
then a suitable form for Q(t) is ¥ A,,¢", where the A,,’s are constants to 
x 
be determined from the boundary conditions. But for any function, 
z = 2(t), of t it is known that if there are two branch points in the region 
with which we are concerned, then in going round one of the branch 
points (across the cut joining the points) the two roots of ¢, as a function 
of z, are interchanged. Thus, in particular, if t is the one root which 
reduces to t, on the boundary, and the other root is e~?*/t then by writing 





Q(t) - Ag+ y A. en a(t" + e-2na/¢n), 
n=1 


i.e. Q(t) = Ay+ ¥ A,(ert"+e-"24-n), (2.1) 


n=1 
it may be ensured that in such a passage round one of the branch points 
the value of Q(t) remains unchanged. 

In effect this is equivalent to saying that in crossing the slit in the z- 
plane from P to P’ say, since these adjacent points P and P’ transform 
to conjugate points ¢, ¢ on the circle |t} = e-*, or tt = e-**, corresponding 
to the slit, then the potential function must have the same value when 
is changed to ¢ where t# = e-®*, 

Thus if O@) = > A, c*", 


x 


we must have 


Ms 


x x = 

% nogn —_ najn —_ —nasg—-n — YS onagn 
> A, 6° = > A, e*8* = 5 A, e-**t-* = 5 A. c*%", 
— @ x — @ 


« 
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so that A_, = A, and the form (2.1) is correct. In considering this, 
however, in terms of the different roots of z = 2(t) it is possible to see 
how the method can be extended to the more complicated transformations 
when there are more than two branch points inside the boundary, as, for 
example, in cases (iii), (iv), (v), (vi), of (IL), section 3. 

This method for determining single-valued potential functions must 
now be applied to the potential functions required for the flexure and 
torsion problem in the particular cases chosen. It is first necessary to 
obtain the boundary values of these functions in terms of t, = e~ rather 
than in terms of { as given in (II). 

The first complete potential Q, is defined (II, p. 7) by 


Q, = (1+-0)Q—oQ = (1+0){b +} —ofbar ta}, 
and in terms of € and » the boundary conditions (I, pp. 86-87) are 


£ fou +80)7],.9 = 0 (2.2) 


and fo = + {o(€)—2(6)}*-+eonst. (2.3) 


(2.2) becomes 


Et =| t ews 20S) ° 
| \ }{2(t) +-2( ya O74 705] 


On ne 


Now "| = |-S] = |" “| = itg| eo] 
7 J ,=0 0€ _ et of _ ot inte 


and thus 
tao ie Ben +2+| 20 x(t) = + ra; (i gL. 


= fe] bot ideo 


so that on the boundary 


2= 2 f fmse()leo-2()} a 24) 


where in (2.4) dashes denote differentiation with respect to t. 
Equation (2.3) may be written in the form 


du = S{an—<()] +-const, 


t=lo 


|& 
~ 


and from this and (2.4) and the fact that 4,, and ¢,, are respectively the 


5002.44 Hh 
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imaginary parts of Q,, and Q,,, we may form our boundary condition, 
viz.: 


Si . l ] 
0,(4,)-0,(7) (1 +2){ Qs (te)—Du()| of Quite) —Ban(7) 


2i[ (1 +0), — ove |-1, 


foo oe) a ao 
: 


(1+e) 
l\ 1 \3 ey 
_of2) — 14 2/4) ef WV sey femal" \lewy_e(! 
Q4 (ty) 0,(7) - ayo) (7) (1 ) | a(t} ( 


== a(t) + 
4 

where dashes indicate differentiation with respect to t. 

Similarly the boundary conditions on Q,, Q, are 


to 
: I Io 3 : fe F = \ 
0,(t,)—0,(7) - Cr ate) +3(7)) i(1+-0) | 208 (7) 2 (t)+42 (7)} 4 


and 2Q.(to) -0,(;] isto(7), 
ty to 


respectively, constants having been omitted. 


“ 


which may be rewritten as 


3. Application of the method to the particular class of cross- 

section 

Let Z z(t) S a, t*, t = et, (3.1) 

n= 

be the transformation which transforms the double Riemann surface 
having a given cut inside the solid section of the cylinder into the ring- 
space between the two concentric circles |f} = 1 and /t e-**, so that 
the two roots of z = 2(t) are t, which reduces to t, = e~*€ on the b~undary, 
and e-?"/t as stated in § 2. 

The constant a, in (3.1) is first adjusted so that the origin is at the 
centroid of the cross-section. The position zg of the centroid is ascertained, 
as in (LI), from 


a | edz ie [ a(7)e0 at 


7 ») Bena Be 
zq8 hi | 22 dz —th 2(03(7)= (i) dt 


where the integrals in the t-plane are taken around the unit circle C,. 
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Bearing in mind the conditions (2.5, 6, 7) we assume that our boundary 
condition takes the form 


l es) 
2(t6)—2(7) = > B,f, (3.3) 
to} n=-—@ 
where we note that B, = —B,. Then writing Q(t) in the form (2.1), viz. 
Q(t) = Ag+ s A, (e"*t" + e-™4-"), 
n=1 


and using (3.3) we find that the arbitrary constants A, must be such that 
A,—Ay = B,, 
A,,e"*—A,e-™* = B 
A,,e-"*—A,, e"* = B_,, 
and thus A,—A,y = By, and for n > 0 


A,, = 4(e"*B, —e-"* B_, )cosech 2na 


n? 


(3.4) 


To determine the coefficients B, in each of the boundary conditions 
we make use of the following notation: 


l @ 

2(t)z|—-}] — ad 3.5 

( (; Stat, (3.5) 

aie l ’ = © > 

2(03(7)= @= > «ef, (3.6) 

n=—@ 

where b= > a,,-4 with t.,= b,, (3.7) 
and c, = > (n+r)a,,,6_,, (3.8) 


which with the origin at the centroid of the cross-section implies that 
cy, = 0. It then follows that 


-(1\.,/1 I< . (l\" 
=(na(7)3 () = — a 2. lj) as ~ 2 &, g-e+8, (3.9) 


3 a 
We also use [2 —a(7)| = > gt, (3.10) 
l 3 wo 
and [:«+a(7)\ = > Af, (3.11) 
n=—@ 
where m= > Sarla,—a,), (3.12) 
and fom F ten ~tedis (3.13) 


r=—@ 
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with a) = a_,. Also h,= > juile_té,) (3.14) 
r= —@ 


and ; 3 b (By 4p +4 _ni7)(@_,+4,). (3.15) 


r= 


4. The three complex potential functions 


The boundary condition for the first complex potential function is, 
from (2.5), 


to 
o,(2) — +24) (NP ey f eayelVeray—e(2\\ 
i) =f se fea 


which with the above notation reduces to 


l ve 
02; (to) -0,(7) > B, nt 
0, r= 0 


42 +} . 
Ww here B, n \ : °) In + \ : es (C,+€ n): 
‘ 12 n 
Thus Q(t) = Ayot > Ayz,,(e™"+e-"%t-") (4.1) 
n=1 


and we have from equations (3.4) 


Aio Ars :(1+20)9,, and forn > 0 





(1-+-20) na —na : 
A,» 7 }(cosech 2n | . ia Jn—e J-n)+ (4.2) 
ot (I a o) f¢ "a(c,+¢_,)+e-"*(c_,+ | 
ee 
Similarly, from (2.6), 
= l «o 
od SB, 
0/ n=—@ 
1 16 
where By» i a ; ed h 2 = teh fe a C_n), 
so that again if . 
Q,(t) = Agot ¥ Az, (e™t"™+e-"t-), (4.3) 
n=1 
we have 
a 19 
Ago—Aa2y = ss +=) ho, and for n > 0 q 
“149 
A, , == }(cosech 2na)| C= (e"*h, —e-"*h_,,)— (4.4) 
“on 
am. — *) fena(c, —&_»)-+e-"%(c C_n —é,)}| 
n 


Similarly from (2.7) 
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so that if Q(t) = Agot > A,,,(e"*t" +e"), (4.5) 
n=1 
then Aso—Ago = tbo, andfor n>0 a 
A;, = hife"*b,, —e-"*b_,, }cosech 2na i 


5. The torsional constant and the centre of flexure 
Following Morris (I, p. 93), we have that the torsional constant I is 
the real part of dQ 
DP 4ar* = bem i) (2G— iz ee, 


but, writing in the notation of (I) 
3 ; dQ, .. 
Z. a Z,—1Z, = u(t, 
we see that r+il* = } | 22(2Z, +72) dz, 
é 
so that 2iT* = 4 [ 23(Z, dz —Z, dz)+}i | 22 d(22), 
é 


c 
and this is zero since on the boundary C the condition 


Z «29 ~ 0 


a ae 
makes Z,dz—Z,dz = 0. 


Thus the integral is real and 


Tl = dp | (292) dz. 
2 z 
This integral being interpreted in the form 


fa@Xy nell) Ln) nal! 
drop | 5 — —i2(7}- coecna(7) dt, 
é 


1 


the genera. esult for T’ is 
r Ampri{ ¥ 2nAy(e"%_,—e-"%,)—4 > Cy Gn}. (5.1) 
= n=—@ 


Following Morris (II, p. 12), but with a slightly different notation, the 
position of the centre of flexure is given by (x9, ¥) where 2, is the real 
part of 


—mpiE-\(A B—H®)-1 ¥ [nb,{AQ, _,—HQ, _,}—nb_,{AQ, ,,—HQ, ,}]— 


n=1 


—|mpE-(A B—H®)-(A—iH)(1—20) ¥ cab, 
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and y, is the real part of 


mpi E-\(AB—H?)+ ¥ [nb,(HQ,_,— BQ, _,)-+nb_,( BQ, , —HQ, ,)]+ 
n=1 


n? 


| Ampi B-(4 B—H?)-\(B+iH)\(1—20) 3 bye 


where £ is Young’s modulus, and Q,,,, Q,_,, are respectively the coeffi- 


cients of t", t-" in the expansion of {Q,(¢)— (1+ 20)z‘(t)}, whilst Q, ,,, 
Q, _,, are respectively the coefficients of t", t-" in the expansion of 
Q(t) — 5 (1+ 20)23(t). 


Also A, B, H are the usual moments and product of inertia of the cross- 
section and are given by 
Pay ne i = 
A+B=-}h | 2z2?dz= —}t | 2(03(7)= (t) dt, 
1 


. . 


Cc Cc 


~ 


. . 


c Cc 


B—A—2H hi | 2*Z dz hi =e(}e"0 at. 


For principal axes, H = 0, and we have that 2, is the real part of 


T7Le( EB) = 2 nb, Qs n nb n! en} = brp( LB)" (1 = 2c) >» b,, C_ n? 
n= n x = 
(5.2) 
Yo is the real part of 


mi(A EB) 2m, n—Nb_,Qy »} + hrpt(AZ)A(1—20) ¥ b,c_,, (5.3) 


n tle x 
and ‘ F 
in | s 44 
A 4° > Cr An g7 x Cy O_n» 
n x n a 
0 «x 
B=jn > ¢,4,—}> > c,2.,. 
n a n= x 
To evaluate Q, ,,, Q,,, we use the notation 
@ 
3 — r = 
” (t) > l, t”, (5.4) 
n XD 
where i + Gitlin 
r=—@ 
x 
and d-= > 6,4,6-» (5.5) 
r= — ao 
and then } A 
Q, n ~~ 1, — (1 -2o)l,; (5.6) 
1 9 ae 
Qs n a A, ,,e**— r(1 T 2a)l,, (5.7) 
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6. The elliptic cross-section 
In the special case of the ellipse we use the transformation 
z = 2(t) = 4cfe%t+e-4-}, (6.1) 
where the double Riemann surface within the ellipse in the z-plane is 
transformed to the ring-space between the two concentric circles in the 
t-plane. 


Comparing (6.1) with (3.1), ie. 2 = y a,,t", we have immediately al! 


the results for the ellipse if we substitute in the formulae of sections 3 and 

5 the values a, = }ce*, a_, = }ce-* 2 
a,=0 (n# +1) 

The required coefficients 6,, ¢,, ete., are 


n? , 
‘ 


by = $¢* cosh 2a, b, = 6_, == }c* (6.3) 
b,=90 (n #0, +2) ' , 
= Ac3e3o, C.,= — Acie 
Cs = 4cke%, Cs = —jcte-* }. (6.4) 
¢,=9 (n # +1, +3) 
f= —9-1 = —Ic(er—e-*)* 
J3 = —J_3 == kc3(e*—e-%)8 e (6.5) 


In =O (n ~ +1, +3) 

h, = h_, = §c8(e*+e-%)8 

hs = h_, = 4c*(e*+e-*)? }. (6.6) 
h,=90 (n~+1, +3) |! 

Since by symmetry the centre of flexure is at the centre of the cross- 
section it is not necessary to evaluate any of the constants for the 
determination of (2, Yo). 

Using these results we determine the constants A,,,, (m = 1,2,3) in 
the potential functions as follows. From (4.2) 


Aj 9—Aio = 0, 
A,, = te{(1+e)(2 cosh 2a+ 1)—(1+-2c)sinh*«], 





4,5 = ~ [ 
~1}8 " 24(2 cosh 2a+ 1) 


A,, =0 (n £0,1,3). 
Q(t) = Ayo+A,(e%t+-e-%t-1) +A, g(e849-+--854-3), 


(1+0)+(1+ 20)sinh*a], 


Thus 


2 8 3c? 
= AyotAj, : 2()+ Aisa 23(t) ——<- ff) 


; 
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and using a = ccosha, 6 = csinha, this result becomes 


ei , i a*-+-a(a?+-6?) , = 
Q(t) = 3a8 p24 (1+-o)+6?}2(t) + 30078) * (t)+real const. (6.7) 
This result agrees with the result given by Love that was mentioned in 


section 2 of this paper. From (4.4) 
As o—Aoo = 9, 
A,, = }ic{(1+ 2c)cosh2a—(1-+-0)(2 cosh 2a—1)], 
ic3 


A, = 1 9 ‘ 2, 2a 1 
- 242 cosh 2a— 1) +20)eosh*a—(1-+0)}, 





A,,=9 (n#0,1,3). 
Thus Q,(t) = Agot+Ao,(e%t+e-%t-1) + Ay g(e® 4? +-e-S4-3), 
which as before reduces to 

—ib? 


. (p21 21 h2)t 
fa?-4267(1+-o)}z(t) ¢ eto +O) 





C006) ase one - nee 29($)-+-real const. 
al?) a?+-3b2' 3(a?+- 36?) (O+ 
(6.8) 
From (4.6), _ 
A39—A39 = Ftc? cosh 2a, Az. = }ic*sech 2a, 
A;, =90 (n # 0,2), 
so that Q(t) = Ag o+As.(e?*t?+-e-2%t-*), 
“(qt?—b2 
or Q,(t) = “7 —"? 22()4 const, (6.9) 


2(a?-+b?)~ 


which agrees with the well known result for the torsion function. From 
(5.1), 


T = hrmpi{4A, .(e?%b_,—e-2%b,)—i(c, 4, +c_, 4_,)}, 
37,3 
that is, fs. (6.10) 
a? +b? 


7. A family of cross-sections with a circular arc as internal cut 

In order to proceed with this cross-section we must search first for a 
transformation of the concentric circle type. The one used here is derived 
from one given by Love (9). 

The transformation is 
4a cot B 


Z z(t) = ~- 
) 2 cot B—i(te*+-t-le-%)’ 





(7.1) 


and the cross-section is one whose internal cut is a circular are of radius a, 


subtending an angle 48 at the point (a,0) in the z-plane, the centre of 
the are E£ being at the point (2a, 0). 
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We may rewrite (7.1) in the form 
c- «Ss 





7.2 
—a kt’ t—ie-*k (25) 


where k = tan}8 < 1 since 0 < 28 < z. 
Remembering now that the cross-section of the cylinder is one of the 
Riemann surfaces only and that we require |t} = 1 on the boundary, we 


z = 2acos pt : 
ve 


ren A ai C 


/ 
i 
a Ze |2a 





-a 





Fic. 1 


choose, as in section 1, the particular Riemann surface corresponding to 
the ring-space e~* < |t} <1 in the ¢-plane. This means that we are 
assuming that « > 0 and that the radius of tie inner circle corresponding 
to the slit is e~-*. Thus, since k < 1 also, we have immediately ke-* < 1 
and the second term in (7.2) can be expanded in powers of (ike~°t-"). 

If further we assume that ke* < 1, i.e. k < e-*, then the first term of 
(7.2) can be expanded in powers of (ike*t) to give as the complete ex- 
pansion 


« 


zs=xsi)= > a,?, 
n=—@ 

where a, = 2acosBi"k"e"*, 
a_, = 2acosPi"k"e-"*, 

The particular form of this cross-section when 8 = }7, e-* = 3, is shown 
in Fig. 1. 

The point C of the figure corresponds to the value t = e~#*” = —i and 
provided the condition e-* > k is satisfied then |z—a| > a at this point; 
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so that OC > 2a. Also OA < 2a under all conditions so that the condi- 
tion e~* > k ensures that the cross-section is the one required with the 
circular are as an internal cut; otherwise part of the cut would lie outside 
the cross-section. 

Next we wish to transfer the origin to the centroid of the cross-section, 
and using (3.2) we find 


2a cos B{k*(e2*— 1)*+-e2%(1 —k?)(1- 1—K*)} (7.3) 


“G 





(1—k4)(1— ke?) (e2@— 2) 

which is real, as was to be expected from symmetry. 
Thus on transferring the origin to the centroid, the transformation 

becomes 





z=2()= > a,t*, 
n ern 


where ad, = 2acosB—z, = 2apcosf, 


~G 


and p —= is a real coefficient which is introduced for 
2a COS | 
convenience; and for n > 0 ( (7.4) 
a,, = 2acosBi"k"e"®, 
a ‘ 2a cos Bi"kne NO } 
All that is now required is the evaluation of the other coefficients b,, c,, 
| n Ch 
ete. From (3.7), 
e2a - 2 
b 4a? cos?B | p? 4 z2I ah as 
° Bit \1—k2e2% ' e2a__ £2) 
and for n > 0 - a (7.5) 
b,, ( )"b .- (—)"d,, 
4a* cos*B|b’en*+ b"(— )"e mee nhn 
where . , Je2e2% l 
7 = Pi —pe 1a’ 
s k2 e2% 
b ss P+ aa Ta e2a 
From (3.8), Ge 0 
and forn >0 c, = 8a° cos*B[c), e"*+-c,(—)"e-"*]i"k”, (7.6) 
where 
, ee n(n—1 nke2™ a 
5°}. sal Eo 
n »2\2 9 22a 62, 2a)2 
(1+?) 2 1—k%e2= * (1 — ke?) 
ys p2a 2 e 
ee G f cacemnmengengay cous equenenaee J aienemene 
(1+ 2)? (1+-k?)? +k ree] * 


+n|p* + k? f= yacte +3: al 


ee <n” ke? e2a 
te (e2x—k2)2 (1+ eta : 
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Also for n > 0 
c_, = 8a* cos®B[c'_,(—)"e"*+-c1_, e-"*]i"k”, (7.7) 
where 
: k2e2% e2m 
c n —_ b 731.©~@632i.9.2 > . > 
(1 —ke23)2 * (1-22 
. } nk* ne e% 
c , _— ———_ es eS 
n (1 | k?)? ’ 1+ k2 1+-e2* (1+-e?%)? 
hn k? —n(n—l) , — ke2« 
ae | de Mee. ek. ANS 
1k 2 TV ta (e22— eau [2 2 
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From (5.4), (5.5), for n > 0 


L, e2na] _ = 24a cos®BT) e™%i"k", 





where 
, 2k? | nk? 2k4 =| (n—1)(n—2) 
- | a, | iss shallhdindeih A lil pated eat ® 7.8 
tn pip+(n If+e) 1+eR 7 +R 6 


The expressions for the potential functions, torsion constant, and centre 
of flexure may now be summarized as follows: 
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Writing, as usual, 2/y = 2(1+<¢) and o’ = a/(1+), we can express 2 
in the form 2 = 2,+0'x 
1 2 
where x, gives the position of the centre of least strain. With the above 
expression for x, we have 
x, B = —2na‘ cos'B(2P—8Q—2R-+-T), 
2, B = —2ma‘cos*p(2P—2R—3T). 

It is obvious that some of the terms in the above infinite series expres- 
sions may be summed to give finite expressions, as has been done in the 
simplest case of Q,(¢). In general, however, finite expressions cannot be 
determined and the rapidity of the convergence of the series depends on 
the values of k and a. For the particular values of k and « used in section 7, 
viz. k = 0-26795, e* = § the series are rapidly convergent, and in parti- 
cular the position of the centre of flexure is given by 

x, =0-1723a, 2x, = 0-0084a 
referred to the centroid as origin. We note that x, is small in comparison 
with x, and the position of this centre of loast strain, as well as the position 
of the centroid in this case, zg = 1-9334a, is marked on Fig. 1. 
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SUMMARY 

The problem of determining the stresses in an infinite strip of finite breadth under 
tension and containing an infinite row of semicircular notches of equal radii placed 
at equal distances on the edges is studied theoretically. Stress concentration is 
calculated by a perturbation method in several cases, the radius of the notches and 
the distance between the centres of neighbouring notches being varied. The decrease 
in stress concentration is studied and compared with the results for an infinite strip 
under tension and containing two semicircular notches placed symmetrically on the 
edges. It is also compared with the stress concentrations in the similar problems 
of an infinite strip or a semi-infinite plate with notches under tension. 


1. Introduction 

[nN 1947 the problem of an infinite strip of finite breadth under tension 
and containing two semicircular notches placed symmetrically on the 
opposite edges was solved by Chih-Bing Ling (1). The most reliable 
determination of the stress concentration in this problem is that of M. 
[sida (2). The knowledge of the stress concentration produced by many 
notches is very important in practical applications, since the stress con- 
centration in this case is smaller than the stress concentration in the case 
of few notches. It is possible, therefore, to improve a design by judiciously 
decreasing the weight of the part. In this paper, the stress concentration 
in an infinite strip under tension and containing an infinite row of semi- 
circular notches of equal radii placed symmetrically at equal distances on 
the edges is obtained and compared with that given by M. Isida (2), in 
each case the radii of the notches and the distances between the notches 
being varied. The stress concentration is also compared with those ob- 
tained by the author (3, 4) for an infinite strip containing two pairs of 
semicircular notches placed symmetrically on the edges and for a semi- 
infinite plate containing an infinite row of semicircular notches under 
tension. 


(Quart. Journ. Mech. and Applied Math., Vol. XI, Pt. 4, 1958] 
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2. Two sets of functions 


Consider an infinite strip of isotropic elastic material under the action 
of a longitudinal tension. Let the strip be bounded in the (x, y)-plane by 
lines equidistant from the z-axis. Let the edges of the strip be notched 
by an infinite number of pairs of equal semicircular notches, with the 
centres of any (arbitrary) pair of notches on the y-axis. Polar coordinates 


T 


suneeeeuee 





(r,@) referred to one of the centres of the notches on the y-axis as pole are 
taken as shown in Fig. 1. The width of the strip, the radius of the notch, 
and the distance between the centres of neighbouring notches are denoted 
by 2b, a (< b), and d respectively. 

Following the method of Ling for solving the problem of an infinite 
strip of finite breadth under tension and containing two semicircular 
notches placed symmetrically on the opposite edges, we will first construct 
two series of biharmonie functions. These series of functions can be 
derived by differentiation from the stress functions for an unnotched strip 
under concentrated (a) normal, or (6) tangential force W acting at each 
of the infinite number of points (0, +6), (+d, +b), (42d, +5),... on the 
straight edges. These two fundamental stress functions are as follows (5, 6): 


2W 
(@) xX%,=—xX 


— 
a 


7 {sinh mb (1-+-mb coth mb)cosh my,—my, 


sinh 5 os ne bait. dm 
\ m — sinh 2mb+ 2mb @  2m* ; 








(1) 
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2W 
Xb = ee CX 
7 
ow x 
* bsinh mb cosh my,—y, cosh mb sinh my, . 
x > mya Ya 44 sin max, dm, 
m/(sinh 2mb-+-2mb) 

is (2) 
where Y= xr+qd, ye=Y, W=7, Yo=y. (3) 


We take as the series of biharmonic functions, which are even in both 
x and y, on account of symmetry, and which give no stresses on the 
straight edges or at infinity, the following: 


Qrr(— 1)R+Ipre+d Grk+2y 





Xx = ; Fy ak+2 
2 W(2k+1)! dx 
inept * (a 
Be 2n(— l yR+1p2k gt, 
Xek ~~ W(2k)! axtett 


where k > 0. To express these functions in terms of polar coordinates, we 
take, in addition to (r,@) already defined, two sets of polar coordinates 
(r,,9,) and (r_,,@_,), where q > 0, and refer to each of the infinite number 
of points (—gqd, —b) and (qd, —b) as poles. Also, for brevity, we express 
the formulae in terms of dimensionless variables as follows: 
a/b=€&, (y+b)/b=y, a/b=A, d/b=8B) (5) 
r/b = p, r, b= Pa § 
Then each function, separated into two parts such that at the poles one 
still possesses the singularities but the other part does not, can be expressed 
as follows: 


any, ' 
Xo = (2k+-1)! | \" +-un)e-¥? + 
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_ 2usinh wy —2u?ne~“ +(1—e~™)(sinh wn —wy cosh un)| - 


‘ 





sinh 2u+-2u 


x u* cos ué(1+-2 s cos qup) du }- (8) 
q=1 


x 
. 





__ 2sinh wy —2une-“— 9(1+e-**)sinh =) * 


, 2 { 
li ioe ania >—UuN) 
Xak aa | \"" sinh 2u-+2u 
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< u* cos ug(1+-2 s cos qup) de | 
q=1 


The relations between the coordinates are 


€ = €+in = p(sin8+i cos 8) 
+ — £+-qB8+in ars p,(sin 0,+-1 cos 8,) (q= +1, +2...) ) 
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and therefore a= $+¢8. (8) 


The first part of the functions x, and 3, can be expressed in closed form 
with the aid of the integrals 


, k! 

| uke-“1 cos ué du = = 008(k+-1)8 
: p 
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x 
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- x Sa 
2 | ute“ cosué ¥ cosquB du - (9) 
qa 1 


0 








— { k! k! 

= > |" cos(k-+ 1)0, + = cos(k+1)6_, 
ont ph P-a@ j 

Further, the results can be expressed in terms of p and @ by using (7) and 


(8). The use of the complex variable greatly simplifies the work and the 
final results are as follows: 





Pa 2k-4 ) eos(2k ee )0, Lp (ak-+ D cos(2k-+ 1)8_, | 
<= ate k+p 2(2k+2p+ 1)! p+? . =" 
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pz 2**® cos(2k-+2)0,+p—2**® cos(2k+2)0_, 
_ S Cpyron__22k+ 2+ Do” 
(2k +1)! (2p)! (gB)**??** 
For the second part the following expansions will be used: 
=: (up)?” 
%n)\! 
hn (2n)! 





cos2p6 (q > 0) 





p=0 


cosh wy cos ug = cos 2n8 


(11) 
. (up)2"*1 

19a | ' 
—, (2n+1)! 





sinh uy cos ug = cos(2n+-1)6 


Hence, omitting trivial terms which produce no effect on the stresses, we 
obtain the following expressions in polar coordinates: 


(2k+-3)cos(2k+1)0 | cos(2k-+3)0 _ 
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By using Dirichlet’s second sida (7), I, and J, may be written 
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uk® uke—4uB 


sinh 4u/B+-4u/B’ fu) = sinh 4u/B+4u/B 





where  f(u) = (17) 





3. The stress function 

When the strip is in a state of generalized plane stress, the stresses 
averaged across the thickness are derivable from a stress function y which 
satisfies the biharmonic equation 


Vty = 0. (18) 
In the absence of the notches, a uniform tension 7' along the strip 
would be given by Xo — }6?T'p%(1-+-cos 26). (19) 


The method of satisfying the conditions when the notches are present 
is to add an infinite set of functions x, and yx, to x, and then satisfy 
the boundary conditions at the rim of one of the notches by adjusting 
the parametric coefficients attached to the functions. Thus the required 
stress function may be constructed in the form 


— (A A 
2 <—2k+1 — eee 
X = xo bP Z (Set xu + pth x) (20) 


Expressing this in polar coordinates, we readily find 
x/b*T = }p?(1+-cos 26) + Q_p?+- 


= 1 [en oe Tae he 
7 Z (4n?—1)]} pe” + 





n=1 
+{(20+ 1)Qon-2 +(22—1)Qan p*}p* |cos 2n0+ 


| S [inte ane Sat IR aa ae 


2n(2n+-1)(2n+ 2) putt 





i) 


R= 


—{(2n—1)(2n-4 2)Qay1+2m(20-+-1)Qan +1 p?}p2"**]e08(2n-+ 1)0, (21) 





where 
~ A A ; 
— 9 fonsa ~ 2k+h 2n+1 an+i 2k 
Qn = 2, (nt Near Ota 

S (2n+1)/i"A, (n > 0) 

cp » (22) 
S » A + n A 

Qon-1 = PS 2n{ (ag, + ane) Se + Base 





— > 2ndz"-1A, (n> 0) ) 
k=0 











484 A. ATSUMI 
and 











2n 2k-+ + 2n+-1 (- 2k- +2n+2 2) Lox 5 on+2— /ok+2n+1—* 2k+2n+1 , 
Ya = ye “92k+2nF1 oe 
2 (—1)*+" \/ 
19 oO]. 
+ (2k+-1) 
. (gp) ten+2 








on 2k+-2n+-2\(2k+-2n+-3)lop.onig |, ‘ 
Pane | 2k+1 Ie a / (2k+-2) 
mn 2k-+2n\(2k-+2n+1)h | 
Ok artis ( 2k iat D2k+2n z+ ant (2k 1) 
52” 1 fase. com 2k 7 2n 1\ {( 2k ‘2 2n-+ -2) Tox zn+2t lop +on +1 1S ok .2n41 
2k+1 sk | D2k+ 2n+l 


~ (—1)* =| / 
oy 2, sored Ihe 
(gB) 2k+2n+2 \ 


A,, are parametric coefficients to be determined from the boundary condi- 
tions at the rim of the notch, that is 


(23) 











o,=0, T= 0 at p=A, —in <I jn. (24) 


The expressions for the stresses are obtained by using (21) and are first 
expanded at p = A into Fourier series between 0 —4tn and 0 = 47 as 
required on the rim. It is found that the resulting expansions contain 
terms of cosines or sines of even multiples of 6 only. Equating the coeffi- 
cient of each term to zero, we have 


: 2QV + fo 0, (25) 
and for n > 1 
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The set of equations may be replaced by 
A gn -gA~™" +- 2NQo ng A*”-* +- (2N— 1) Qo, A" + (2n—1)F,, +48, = 9, 








(28) 
G,, + Von-2 pe. 2Qon A" + Qon sg A*"t2+ 48) 5 = 
where eae 
F,, Dd ("be5—2 A op-1A-* 1 — "pg Dog A*-*) 
s=1 > (30) 
G, = p> ("Ypog—2 Aggy A“? — "Vg Dag A*-*) 
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a2 n{4s?—(2n—1)*}(2n+ 28+ 1)’ wt inee—1 é 
eh be 8(—1)"*+4(2s—1) See 
Mes {4s*—(2n—1)\(2n—2s-+ 1)’ * an—e—1 ™ 
(31) 


These equations for the biharmonic stress function and boundary condi- 
tions remain invariant even if the equations are expressed in terms of 
polar coordinates associated with any other one of the notches. Then, 
if the boundary conditions on one of the notches are satisfied, the corre- 
sponding conditions on the other notches are automatically satisfied. 

A formal solution of equations (28) and (29) can be obtained as in 
Ling’s paper (1). But, since this method is one >f laborious successive 
approximations, a perturbation method, as in Isida’s paper (2), in which 
A is the perturbation parameter, is introduced here. 


4. A perturbation method 


In order to expand A,, into series of Q,, we first require to solve the 
sets of equations for the (A,,,,,A~*"-*) given by (29) after setting Q,A* = 1, 
Q,,A" = 0 (m #£k). Then substituting the given values of (A,,,,,A~*"~*) 
into (28), the values of (A,,,A~?"-*) are to be obtained. These calculations 
are not carried out here since Isida (2) gave the results, valid for this case, 
by solving the eight sets of equations obtained by setting k = 1 to 8. 
When the values of (A,,A~"-*) are thus obtained, the values of (A,,A~"-*) 
can be determined as linear multiples of Q,A*. The relations between A,, 
and Q, are then expressed in the form 


Ay, = (QA *)aR-+ Y (—1)A#3faB/(k-+1)} Op, (32) 


where the values of a? are the same as shown in (2), p. 8, Table 1. Next, 
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assuming that A, and Q,, can be expanded in powers of A when A is small 
the expansions may be obtained from (22) and (32) as follows 


A - Ss Aln+pt+2)\n+p+2 
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Equating the terms of the same order of A in both sides after substituting 
(33) into (22) and (32), we obtain 
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0 at the rims of the notches, the maximum stress is given in 
the form 
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Substituting (33) into (37) determined by using (21) and (22), and then 
expanding in powers of A, we obtain 
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(n+6 (3), gp(2)\y4 
p 3 Ar — Oi" + O3 jA T 
n=0 
LIpow_ +7 (4 (2)_1 (3)! 5 
\Vo > Ay —Q} —Q§ +. O§ jA° 
From (36) and (38) the numerical tail for o 


+... (38) 
(8). . 


T can be obtained for 
any given value of 8. When A = 0, the results are equivalent to Maunsell 





J 8 
When A + 0, it is necessary to prove convergence of the series 
though this is usually difficult 


However, by comparison with Howland’s 

















STRESS CONCENTRATIONS IN A STRIP UNDER TENSION = 487 
work (9) and from the agreement of the numerical results with experi- 


ment, convergence looks likely for A < min(0-5, 0-258). The numerical 
calculations are made in the lower parts of this range. 


5. Numerical calculations and results 

The numerical calculations have been carried out for the cases 8 = 0-5, 
1, $7, and 2. First an evaluation of the integrals J, and J, was performed 
in each case. Then the values of y3?, y3%.,, 53271, and 8327} were deter- 
mined by substituting the values of J, and J, into (23). Finally, o,,,/T 
was obtained by using (36) and (38). In each numerical equation, for 


2 «x 
the terms of A, A", and Al, — ¥ az Q!/9, — > (a3 QL/9—aZ Qf?/10), 
n=0 n=0 
@ 
and — > (az QY/9—az QY/10+-a}, O{B/i1) were omitted respectively. 


But, olan the terms omitted seem to be practically negligible, the following 
equations may give the results to an adequate degree of approximation. 
They are 
Omax/ 1(B = 0-5) = 3-065 —(9-7330 x 10)A?+-(1-3781 x 10-5)A3+- 
-+-(5°9708 x 108)A4+-(8-3745 x 10-4)A5— (2-8128 x 105)AS+- 
-+-(2°2097 x 10-2)A7 +-(1-2386 x 107)A8+-(3-7545 x 10-7 )A®*— 
— (52588 x 108)A!+-(3-3195 x 10)A™ +-(2-1880 x 101)Al, 
(39) 
Omax/ L'(B = 1) = 3-065—(1-6733 x 10)A?+-(2-4696 x 10-*)A3 +- 
|. (2-6955 x 102)A4+-(4-8029)A5— (2-7816  103)A8+- 
+ (3-6025 x 1O)A7 + (2-6392 x 104)A8+-(1-2558 x 102)A9— 
— (2-7031 x 105)A!9+ (6-6733 x 102)Al-+ (28017 x 108)A12, 
(40) 
max! 1'(B = 47) = 3-065—(3-9658)A2-+ (3-8920)A9 + (1-6983 x LO)AS+ 
4. (37174 x 1O)A5— (17500 x 10?)A®+-(9-0669 x 10)A7 + 
+. (2-8368 x 102)\8—(1-5480 x 10)A9—(9-3469 x 102)A10— 
| —(1-8312 x 103)A™ + (8-9066 x 103)A?2, (41) 
‘ii T(B = 2) = 3-065—(2-2238)d2+-(8-2007)A3—(1-6238 x 10)A*++ 
4.(5-1573 x 10)A5—(7-8374 x 1O)AS+ (2-3089 x 10)A7 + 
4.(1-5464 x 102)A8—(4-5787 x 102)A®-+- (1/1430 x 103)Al0— 
— (24943 x 103)A1 +4 (3-7532 x 103)A12, (42) 
In Table 1, the values of o,,,,/7' calculated from the above equations are 
shown and are compared with those given by Isida (2) for Ling’s problem. 


/ 
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TABLE 1 
Values of O/T 
A B=o'5 B=1 B=42 B= 2 M. Isida 

oo 3°065 3°065 3°065 3°065 3°065 
0°05 2°86 3°03 3°06 3°06 

0075 2°67 2°98 3°04 3°06 

o'lo 2°92 3°03 3°05 3°05 
ors 2°80 3°00 3°04 

0°20 2°70 2°07 3°03 3°94 
or25 2°95 703 3°04 
0°30 2°95 3°03 3°95 
0°35 3°06 

040 311 

3-0 
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In Fig. 2, the values of stress concentration factor x = o,,,,(1—A)/T 
calculated from the values in Table 1 are given. Here those given by 
Isida are shown in dotted line. In Fig. 3, the values of o,,,,/7' in Table 1 
are compared with those calculated by the author (4) for a semi-infinite 
plate under tension and containing an infinite row of semicircular notches. 
The latter is shown in dotted line. In Fig. 4, the values of the stress 
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concentration factor x = o,,,(1—A)/7' given in this paper for 8 = 1-0, 
those given by the author (3) for an infinite strip in the same conditions 
but only containing two pairs of semicircular notches, and those of 


M. Isida (2) are compared. The values of the last two are shown in dotted 
lines. 


In conclusion, I wish to record my indebtedness to Professor Emeritus 
S. Higuchi of Tohoku University for much kind advice and many stimu- 
lating discussions. 
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THE FLEXURAL VIBRATIONS OF A CUT THIN RING 


By L. 8. D. MORLEY 
(Royal Aircraft Establishment, Farnborough) 


[Received 3 October 1957] 


SUMMARY 


There appears to be little published work on the flexural vibrations of a cut thin 
ring apart from a recent paper by H. Hasselgruber (1) who used the Rayleigh—Ritz 
method to find the frequency and shape of the first three modes. In the present 
paper these vibrations are investigated using exact methods, and tabulations are 
given for the first ten modes of symmetrical and anti-symmetrical vibration. There 
is no difficulty in extending the method to treat many problems of the flexural 
vibrations of circular ares. 


Nomenclature 
We shall use the following notation: 


T = circumferential tension. 

Q = radial shearing force. 

M = bending moment. 

R = radius of the ring. 

A = cross-sectional area of the ring. 

w, v = radial and tangential displacements respectively. 
m = Mass per unit circumferential length of the ring. 
k = radius of gyration of the cross-section of the ring. 

t = time. 
@ = angular coordinate, origin diametrically opposite the cut. 
E = Young’s modulus. 


i 


1. Fundamental equations 


THE coordinate system for the cut ring is shown in Fig. 1, and the forces 


acting on an elemental portion are shown in Fig. 2. For small displace- 
ments the equations of motion are 





aQ ? ew ) 

a T= mR—s 

eT ov 

a9 + 9 = mR ee i 
eM eer ow 
“ero 5”— a0) 


(Quart. Journ. Mech. and Applied Math , Vol, XI, Pt. 4, 1958] 
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oo, 
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etl 
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Fic. 2 


the relationships between the strain and tension and bending moment 


7 — FA 4) 


being respectively 


R \e6e 
(2) 
EAk* c ew 
M = 
R = G{e- 3) 


Now, it is usual in this type of analysis to neglect the rotatory inertia 
effect, which amounts to neglecting the right-hand member of the third 
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of equations (1), and to assume that the central line of the ring is inex- 
tensible. From the first of equations (2) this implies that 


ov 
= ati aaa 3 
w 56 (3) 


Resulting from these considerations, the governing equation for the motion 


is found to be 
(F/,, CY _ oh, et 4 
| + 5m otk? at w)| gin. (4) 


where c is the velocity of wave propagation in a thin straight rod, i.e. 
| EA 
c= j—.- 
m 


The normal functions for free vibration in the symmetrical modes are 
determined by taking the complete primitive of equation (4) to be of the 
form 


2. Symmetrical vibrations 


3 
M = ¥ a,.cos(A,, 6)cos( pt +-e) (5) 
xk=1 


where the a, ar, as yet, arbitrary constants, p/(27) is the frequency of 
vibration, and the A, are the roots of the equation 


Rp? _ d(1—D*)? 





- — > . 6 

c*k? 1+-A? : " 

The corresponding expressions for the circumferential tension and radial 
shear are 3 

_ zy a x coa(A @)cos( put +-€) (7) 

~ 3ST 

\< 

Bas — BD, %-Acsin(A, Soos(pt-+), (8) 


and for a complete record, the radial and tz ngential displacements are 


3 


ee: R Ber 
u — FAR PS a ai @)cos( pt+-e) (9) 
, R? S Sere 
anc v= BAB 2, “yi )eos( pt +-e). (10) 


The frequency p/(27) becomes determinate when the boundary condi- 
tions are satisfied and this requires the vanishing of the bending moment, 








494 L. 8S. D. MORLEY 


circumferential tension, and radial shear at @ = 7. This gives the deter- 
minantal equation 


l l l 
/A(AZ—L)jttand,7 A,(AZ—1)?tanA,7 A,(AZR—1)2tandAg7 |= 0 (11) 
| @j-1)" (3-1)? (A3—1)? 

which is valid provided that there are no equal roots A. 

In the evaluation of the frequencies from equations (6) and (11), it is 
soon found that there are advantages in using A, as the independent 
variable, so that A,, A;, and p are dependent variables. Now, equation (6) 
is a cubic in A? and we may say that Aj must always be real and, further- 
more, lies outside the range —1 < Aj < 0 for real and non-trivial p. 


Again, from this equation, it is readily deduced that 








anni err R'p? 

A? \3+-AZAZ+-AZ A? — 1 - cape | (12) 

R'p? 

222 — 
MAGS = is 
+ xe. 2-0, 4(1—A2)?>)h 
so that Ag, Ag = ——2 14 {[1-——__._+,, ;}" |, 13 
| £8 ="S]('- T pea ity 
and from this last it is seen that if A? < —1, then either Aj or A2 is positive. 
Hence, there is no loss in generality by considering only values 

Aj > 0, (14) 


and it is worth noting that when 0 < A? < 3(7+-~17), the roots A3 and A3 
are complex and an approximation to A, gives a purely imaginary value 
to the determinant of equation (11). 

Expansion of equation (13) yields 


= M+3+5¢~ | 


4 
P mh (15) 
Mi 
and so the determinantal equation (11) reduces to 
A, tandA,7 = A, tandA,7 (16) 


when terms of order A; * are neglected in comparison with unity. Further- 
more, substituting the first of equations (15) into this last, it is found that 


3 
tanA, 7 = +o 


3 
or, in other words, A, = (n—}) 4+ 
1 Ost Ge 
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where n is an integer. This formula improves in accuracy as the integer n 
increases. The values of A,, Ay, As, and R*p/ck have been calculated for 
the first ten modes and are appended in the table below. 








TABLE 1 
Mode A,t Ay As R*p|ck 
I 1°24927 0°53387+0°255591 | 0°53387—0°255501 "043767 
2 1°85206 0°46878+ 0966861 | 0-46878—0-966861 2°13835 
3 2°78662 2°15403t 1-060841 6°36764 
4 3°76855 3°343078 1°01285¢ 12°7604 
5 4°76122 4°433971 1°00452t 21°2065 
6 5°75752 5*490451 1°002011 31°6749 
7 675540 6529421 1°00103t 44°1543 
8 7°75407 7°558071 1-000587 58-6399 
9 875318 8: 580061 1°000351 75°1294 
10 9°75255 9°59749 1090231 93°6213 

















+ These values of A, were calculated by Dr. K. N. Dodd, Royal Aircraft Establishment. 
t Hasselgruber (1) calculated these values as 0-4377 and 2-1856 using the Rayleigh- 


Ritz method. 


The values of the arbitrary constants a, in equation (5) are given by 
a, = B(AZ—A3)(1+-Aj)sec A, 7 
a, = B(A3—A?)(1+-A3)secA, 7 }, (18) 
a, = B(A{—Az)(1+-Aj)sec Ag 7 

where § is a normalizing constant such that 


[ M?d0 = cos*(pt+e). (19) 
It is worth noting that since the bending is the only contribution to the 
potential energy then 


[ M, M,,d0 = 0 (20) 


for any mode n not identical with mode m. The calculated values of the 
normalized constants a, are given in the table below. 








TABLE 2 
Mode a ay as 
I 02612 0°2019+ 0°13821 o°2109—0°1 3822 
2 —o'5256 0°004233+0°022271 | 0°004233—0°022271 
3 0°5557 01053 X 1078 —0o'1568 x 10-% 
4 —o5612 2313 X107* —0°1477 X 107% 
5 05628 07352 x 107% —0o°3041 x 10~¢ 
6 —0'5634 0°2633 x 1077 —0°8934 x 1075 
7 0° 5637 o 1001 x 10-8 —0°3264 x 107% 
8 —0'5639 03937 x 107!” — 01386 x 1075 
9 05640 01584 x 10-4 — 06565 x 107% 
10 —o5640 06467 x 1078 —0°3384 x 107* 
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3. Antisymmetrical vibrations 

The analysis for the antisymmetrical vibrations follows precisely the 
same lines as just discussed and so a catalogue of results will suffice. 

The bending moment, circumferential tension, and radial shear are given 
by the expressions 


3 
M = > a, sin(A,,@)cos(pt+-e), (21) 
K=1 
2 2 
T —*% 2% = 2 ne #)cos( pt+-e), (22) 
I~ 
Y= BZ, 8-0. Poe eh, (23) 


and the radial and tangential displacements by 


w par 2," : ra ind, #)cos( pt-+-e), (24) 
R2 Z, l 
e EA k2 Ps Ge Las) Oe Poos( Pt €)- (25) 


The determinantal equation for the evaluation of the frequencies p/(27) is 
l l l 


A,(AZ—1)? cot A, 7 A,(AZ—1)2 cot A, 7 Ag(AS—1)? cot A, z 
(\2—1)2 (AZ—1)2 (A2—1)2 


| 
i= Q (26) 


which is valid provided that there are no equal roots A. The same remarks 
apply for the evaluation of the frequencies except that the determinantal 


equation (26) reduces now to 
A, cotA, 7 = A, cota, 7 (27) 


, 
when terms of order A;* are neglected in comparison with unity. The 
approximation for the root A, is now 
3 


—_____— 28 
4(n—})?’ = 


A, = (n+}4)+ 


where » is again an integer. This formula improves in accuracy as the 
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integer n increases. The values of A,, A,, A; and R*p/ck have been calculated 
for the ‘first ten modes and are appended in the table below. 








TABLE 3 

Mode A,t Ay As R*p/ck 
I 146556} | 0°53527+0°600371 | 0°53527—0°600371 0°94816§ 
2 2°30773 o 16891 + 1°30052% 0° 16891 — 1°30052% 3°96902 
3 3°27532 2770671 1°025221% 9°30373 
4 4°26419 3°89 4691 1°007 33% 16°7295 
5 525909 4965091 1-002951 26-1888 
6 6°25633 6-011 56% 1°001 421 37°6636 
7 7725466 7°04476t 1°00077t 51°1465 
8 8-25358 8069741 1'00045t 66°6343 
9 9°25284 9°08g25% 1000281 841251 

10 10°25230 10°104921 1'OOO19t 103°618 

















+ These values of A, were calculated by Dr. K. N. Dodd, Royal Aircraft Establishment. 
t This value is accurate only to within +0-00002 and A,, A, and R*p/ck have been 
calculated from this value. 
§ Hasselgruber (1) calculated this value as 0-9709. 
rhe arbitrary constants a, in equations (21) to (23) are given by 
a, = (A}—A§)(1+-A})cosec A, z, 
(A3—Aj)(1+-Aj)cosec A, 7, 
a, = (AZ—A§)(1+-A})cosec A, 7. (29) 
It has, unfortunately, not been possible to tabulate the values of the 
normalized a, for these antisymmetrical vibrations. 


l 


a, 
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TRANSVERSE WAVES IN ELASTIC PLATES 


By V. T. BUCHWALD 
(Manchester College of Science and Technology) 


[Received 26 September 1957] 


SUMMARY 

Expressions are derived for the phase and group velocities of plane rotational 
waves in elastic plates with displacement vector parallei to the surfaces of the plates 
and perpendicular to the direction of propagation. It is shown that in a double 
plate two types of transverse wave can travel along the plate, each type consisting 
of an infinite number of modes. The behaviour of the waves is illustrated by 
dispersion curves, computed for the lowest modes. The paper concludes with an 
investigation of the properties of Love waves in a surface layer, including a 
discussion of the behaviour of the minimum group velocity, discovered by Jeffreys 
(1), for which a lower bound is found. 


1. Introduction 

ReEcENT work on the propagation of vibrations in elastic plates has been 
almost confined to the study of displacements in planes parallel to the 
direction of propagation and at right angles to the faces of the plates. 
Whilst it is true that most of the energy of a disturbance will be trans- 
mitted by waves of this form, a proportion will, in general, travel as waves 


of rotation about axes perpendicular to the faces of the plates. Some pro- 


perties of waves of this second kind were investigated by Love (2), Jeffreys 
(3), and Stoneley (4), in connexion with the propagation of waves on the 
surface of the earth, but their presence in plates has been neglected so far. 

If it is assumed that at a large distance from a disturbance waves of 
this kind travelling along a plate have plane fronts, their displacements 
will be parallel to the faces of the plates and transverse to the direction 
of propagation. It will also be assumed that these waves can be resolved 
into harmonic components so that with the use of harmonic analysis 
results of a very general nature are obtained which are applicable to most 
kinds of initial value problems. 


2. The equations of motion and the boundary conditions 

A rectangular cartesian system of coordinates O(x,y,z) is chosen so 
that the faces of the plates are planes z = const, and (u,v,w) denote 
the x, y, z components of the displacement vector u. In the case of a plane 
wave travelling in the direction of the x-axis the z-component of the 
rotation VAu, a 


Cx Oy 


(Quart. Journ. Mech. and Applied Math., Vol. XI, Pt. 4, 1958] 
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will correspond to a displacement in the plane of the faces of the plate 
and perpendicular to the direction of propagation. The equation of 


motion is a2 
Ory 
272), — 
civ = —, 
, ot? 


where c, = (u/p)! is the velocity of rotational waves, p is the density, and 


(2.2) 


» is Lamé’s constant. 
The conditions to be satisfied on a stress-free face are 


~~ 


22 = 22 = zy = 0. 


From the stress-strain relations we find that the condition on y¥ at a stress- 


free face z const is ~ 
op _ 9 


“= (2.3) 
cz 

If the plane z = const is the boundary between two media with different 
elastic properties then the stresses and displacements must be continuous 
across that face. This means that % and pu(é/éz) are continuous across 
the boundary. 

For the remainder of this paper we confine our attention to plane waves, 
v being the only component of the displacement considered. 


3. Transverse waves in a plate 
The plate is assumed to occupy the space between the planes z = +h 
with boundary conditions 
yo on z = +A. (3.1) 
oz 
Let % be composed of harmonic components of the form 
yb ‘ain Weilpr =n) | (3.2) 
representing a plane wave travelling in the direction of the x-axis. V is a 
function of z only, n/27 is the frequency, and A = 27/p is the wavelength 
of the waves. Substituting for y% in the equation (2.2) gives 
GY (n* - 
EM a inn (3.3) 
dz* ' \@ 
where c, is the velocity of rotational waves in the medium. The general 
solution of the differential equation (3.3) is 
VY = Acosrz+ Bsinrz, (3.4) 
where r? = n?/c?—p?. (3.5) 
Thus, using the boundary conditions (3.1) it is seen that the only possi- 
bilities are that either A = 0, rh = (m+4)z, or B= 0 and rh = mz, 
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where m is any positive integer. If A = 0 the displacements will be anti- 
symmetrical about the plane z = 0, while in the case B = 0 the displace- 
ments will be symmetrical. We shall devote our attention to the anti- 
symmetrical case, results for the symmetrical case being entirely analogous. 
For antisymmetrical displacements 


ab — Beipr- nt) gin Fath (3.6) 


where Tm = (m+4)r/h (3.7) 


m 


and n®/ci—p* = r%,. (3.8) 


Thus ¢ consists of an infinite number of branches, or modes, depending 


first mode, that with m = 1 the second, and so on. From the equation 


(3.8) it is clear that n? > (7»)?, (3.9) 


so that the minimum, or ‘cut off’, frequency in the (m+ 1)th mode is given 


by n = (m-+-4)me,/h. (3.10) 


Thus we see that no wave of this type can be propagated for which 
nm < ¢,/2h. 

Writing c (= n/p) for the phase velocity, we see from (3.8) that 

¢ = ¢,(1+15,/p*)}. (3.11) 

As the ratio p/r,, varies from zero to infinity, v varies from infinity to the 
limit c, for waves which are very short compared with the thickness of the 
plate. 

The group velocity V (= dn/dp) is found from (3.8) to be given by 


V = e,(1+12,/p?)-*. (3.12) 


Thus for short waves p/r,, is large and V will tend to the upper limit c,, 
while for long waves p/r,, is small and V tends to zero. 


It is beyond the scope of this paper to discuss at length the concept 
of group velocity and good accounts may be found in Jeffreys and Jeffreys 
(5) or Stratton (6). The disturbance associated with a particular wave- 
length will travel with the group velocity, while individual waves will 
travel with the phase velocity, but will change their periods, lengths, and 
velocities as they travel. In view of the simplicity of the expression for 
the group velocity in (3.12) it is possible to discover what happens to a 
wave after a long time. The total disturbance travelling as transverse 
waves may be expressed in the form 

o(2, 2,¢) = » sinr,,, 2 [ A,,(p)et2-™ dp, (3.13) 
m= ¢ 


x 
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where A,,(p) depends only on the initial conditions at t = 0. Sufficient 
conditions for v and & to exist in this integral form are that 


A,,(p) = O(1/mp*), (3.14) 
as m or p tends to infinity. The integral in (3.13) can be readily evaluated 


for large x and t by the method of stationary phase, which yields the 
asymptotic expression 


v(x,2z,t) ~ (2m)teyt > rt sinr 2A lest, (3.15) 

éi = é 
where £ = c?(2—x*. Remembering to take the correct branch of é+, it 
may be seen from the condition (3.14) that when x = c,t, v is zero, while 
if x > c,t, p becomes imaginary and the exponential term in (3.15) tends 
rapidly to zero. Thus it follows that the maximum velocity of propaga- 
tion of the disturbance is c,, and that the shorter waves will arrive shortly 
after the first signs of the disturbance, while the longer waves will suffer 
considerable dispersion and will travel more slowly. In addition, it can 
be shown that for a fixed wavelength the lower modes will suffer less 
dispersion and will travel more rapidly than the modes for which m is 
large. 


4. Transverse waves in a double plate 

We shall assume that the plate consists of two different materials, one, 
in which the velocity of rotational waves is c,, occupying the space 
0 <z<h, while the other, in which the velocity of rotational waves 
is C,, occupies the space —h < z < 0. There is no loss of generality in 
assuming that c, is greater than c,. Here pu, w, are the respective Lamé’s 
constants for the two materials, and 7 is the ratio p/p. 

If %, and &, are the respective displacement functions in the two layers 
then the equations of motion are 

Oy, Oy _ 1, 


atta awe when 0 <z < Ah, (4.1) 


and O's Obs = kt Orbe when —h <z < 0. (4.2) 
dat ' Os Gk OF 


The boundary conditions are 


 . onz =h, 
z 

ss on z = —h, 
ez 


Uy = Yay wy(2py/@2) = pg(2pg/@2) on z = 0. 
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We assume that f, = Y, efws-n0, (4.6) 
and that §, = Y, cor-™, (4.7) 


Substituting ¥, and ys, from (4.6) and (4.7) in the equations of motion 
gives the equations _— : 
© d? 1, n= 

atts WE 2p — 0. 4.8 
dz* (" : , 


d?¥, n? 
and +t 
dz? 


If c = n/p is the phase velocity then there are three possibilities. Either 


ae 0. (4.9) 


9 
C3 


(i) ¢ >), or (ii) c, > ¢ > Cg, or (iii) c, > ¢. 

Stoneley (7) treated the problem of a double plate as a limiting case 
of a double surface layer, but did not consider the possibility c, > c. He 
therefore erroneously concluded that the wave velocity in a double plate 
is greater than the velocity of rotational waves in either of the two layers. 

(i) Taking c > c, we put 

r? = n*/c?}—p?, and rs = n®/c3—p*, (4.10) 
where r, and r, are real positive numbers. In this case the general 
solutions of the differential equations (4.8) and (4.9) are 

", = Acosr,z+ Bsinr,z, (4.11) 
and Y, = Ceosr,z+Dsinr,2z, (4.12) 
where A, B, C, and D are constants. If these expressions are substituted 
for ¥, and , in the boundary conditions (4.3) to (4.5) and the arbitrary 
constants are eliminated we obtain the transcendental equation 


mr, tanr,h = —r,tanr,h. (4.13) 


This may be taken as an equation connecting the dependent variable n 
and the independent variable p, with the solution for » consisting of an 
infinite number of branches corresponding to modes of vibration. Since 
{ 


1 > €, it follows from (4.10) that r, > 7,. If, lies in the interval 


0 < hr, : hr, 


then it is evident from (4.13) that r, is at least }h7. Thus we see from 
(4.10) that the minimum frequency of this type of wave is given by 
n > ¢,/2h. 
The phase velocity c is obtained from the equations (4.10), which yield 
ce? = ye}(r}—r})/(r3—yr}), (4.14) 
where y = cj/c3. Differentiating the equations (4.10) and (4.13) with 
respect to p we obtain the group velocity V in the form 


V = ye}(or,+r,)/e(ydr,+r2), (4.15) 
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where th = (hr, sec*hr,+-tan hr,)/(hr, sec*hr, + tan hr,). (4.16) 


The shape of the dispersion curves of the phase and group velocities will 
depend on the branch of the solution of (4.13) that is taken, as well as 
the values of the constants y, r. Although no solution of these equations 
exists in closed form the equations (4.10) and (4.13) can be solved 
numerically to give corresponding values of r,, r., p, and n, whence c and V 
can be evaluated from the equations (4.14) and (4.15). If we devote our 
attention to the lowest mode in which 0 < r,h < 47 and $7 < r,h <7, 
we see that as r, tends to zero, hr, tends to 7, c to (y)!c,, and dr, to m/2hr, 
so that V tends to (y)#(1+-2r)e,/(y+2r). In addition, ifr, = 0, 
(hp)® = #*/(y—1), 

and this expression gives the maximum value of p in this mode. As r, 
increases we see from (4.14) that the largest possible value of r, occurs 
when r2 = yr?, in which case p = 0, ¢ = 0, and V = 0. These features 
are illustrated in Fig. 1, where computed values of c and V are plotted 
against values of ph, the constants y, + being given the values of 1}, 1 
respectively. The most interesting feature of the graph is that the group 
velocity is found to have a maximum of 1-06c,, when ph is approximately 
4. Jeffreys [(5), pp. 515-18] has shown that when the group velocity 
has a stationary value, the dispersion will be comparatively small and 


that comparatively large amplitudes will be found at a large distance 
from the source. Bearing this in mind, in this particular case we would 
expect the front of the wave to travel with a velocity of 1-06c,, with wave- 
lengths near the value given by ph = 4. The longer waves would arrive 
later and would have comparatively small amplitudes. 


(ii) Case c, < ¢ <¢y. In this case we take 

r? = p*—n?/c?, r3 = n*/c§—p?, (4.17) 
r, and r, being real and positive. Here the general solutions of the 

equations (4.8) and (4.9) are 
Y, = Acoshr,z+ Bsinhr, z, (4.18) 
and Y, = Ceosr,z+ Dsinr,z. (4.19) 
If these equations are substituted in the boundary conditions we obtain 
the equation rr, tanhr,h = r,tanr,h, (4.20) 


which again is an equation for n in terms of p, the solution possessing an 
infinite number of branches. As the left-hand side of (4.20) is always 
positive, r, in the mth mode must lie in the interval 


(m—1)2 < hry < (m—4)n. 
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= Phase velocity 
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From the equations (4.17) we see that as Ar, tends to (m—1)z, r, tends to 
zero, so that p tends to a lower limit given by p? = (m—1)*x?/(y—1)h? 
while the minimum ‘cut-off’ frequency is given by 


n? = y(m—1)*x%c,/(y—1)h?, (4.21) 


for the mth mode of vibration. When hr, tends to the upper limit (m—4)z, 





TRANSVERSE WAVES IN ELASTIC PLATES 505 
both n and p become very large. From (4.17) the phase velocity c can 
be expressed in the form 

f= yet rebl(t +r). (4.22) 
The group velocity V is again given by the expression (4.15), but in this 
— th = (tan hr,+hr, sec‘ hr,)/(tanh hr,+Ar, sech*hr,). (4.23) 


In all modes except the first we see that if hr, = (m—1)z then r, = 0, 
c = yes, ror, = 4r, 80 that V = (1+ 27)y*c./(y+2r), whilst if hr, = (m—}4)z, 
r, is infinite, so that both c and V are equal to c,. Thus all the modes 
above the first will have similar properties, which will depend on the 
values of the constants y, 7. The first mode is exceptional because, from 
(4.20), r, = rr, when r, is very small. From this it can be shown that 
as r, tends to zero the limits of c and V are both [y(1+-7?)/(y+7*)]te., 
while, as hr, tends to the upper limit 47, c and V tend to the limit c,. 
The properties of the first mode are illustrated in Fig. 2 by graphs of c 
and V against ph, where, for the purpose of computation, the values of 
y, tT have again been taken as 1-5, 1. The main features are the maximum 
at p = 0, when V = c = 1-095c,, and the minimum group velocity around 
ph = 5, where V = 0-989c,. 

(iii) Case c <c,. Here the solutions for ‘Y, and ‘Y, are hyperbolic 
functions, which lead to the equation 


tr, tanhr,h = —r,tanhr,h, (4.24) 


when combined with the boundary conditions. A glance at this equation 
shows that it cannot be satisfied by real values of r, and r,. This means 
that no wave can travel along the plate with a phase velocity less than cy. 

Thus transverse waves of two kinds can be transmitted along a double 
plate. If the thickness 2h of the plate is large compared with the lengths 
of the waves, then it can be seen from the equations (4.6), (4.7), (4.11), 
(4.12), (4.18), and (4.19) that the waves of the first kind would correspond 
to transverse waves travelling through the material and being refracted 
and reflected at the face z = 0, as well as reflected at the faces z = +A. 
The waves of the second kind now correspond to transverse waves 
travelling in the material z > 0, being reflected back at z = h and ‘totally 
reflected’ at z = 0, with the vibration only penetrating a very short 
distance into the other medium. 


5. Love waves in a surface layer 

We consider a surface layer 0 < z < h, in which the velocity of rota- 
tional waves is cy, resting on a semi-infinite material which occupies the 
space z < 0 and in which the velocity of rotational waves is ¢,. 
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Here the equations of motion are the same as in (4.1) and (4.2), and we 
take ys, and ys, in the form 


yf, = (Acosr,z+ Bsinr,zje?-™, mo<z<h, (5.1) 
and ys, (Clenzt+i(pr nt) in z < 0, 


where r, and r, satisfy the equations 


ri = p*—n/c? and r3 = n?/c3—p?, 


so that n® = yc3(r? +-r3)/(y—1) ) 


and p? = (yr? +r2)/(y—-1) I 


where y = c}/c3. 
The boundary conditions are é,/éz = 0 on z = h, while %, = 4, and 
7 Oy, /€z = Oy,/ez on z 0. Substituting the expressions for %, and ys, 
from (5.1) and (5.2) in the boundary conditions gives the equation 
7, = r,tanr,h, (5.5) 
first derived by Love (2). This may be regarded as an equation for n in 
terms of p with the solution possessing an infinite number of branches. 
In the mth mode r, lies in the interval given by the expression 
(m—1)z < hr. < (m—})z, 
and a glance at (5.4) shows that in this mode 
(hp)? > (m—1)*x?/(y—1), 


and (hn)? > ye3(m—1)?x?/(y—1), 
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the latter inequality yielding the cut-off frequency for each mode. The 
phase velocity c may be obtained from the equations (5.4) in the form 


ce? = yes(ri+r3)/(yri +13), (5.8) 

and if we differentiate (5.3) and (5.5) we obtain the group velocity 

_ [erin thr (+77) yes 

ore tryin it epe 
If we look at the mth mode we see that if hr, is at the lower limit (m—1)z 
then r, = 0, and both ¢ and V tend to the upper limit c,. In particular, 
we see from (5.6) that the longest waves only occur in the first mode, in 
which we see from (5.5) that rr, = Ar? as r, tends to zero. Using this 
result in (5.8) and (5.9) we see that for small r, 





(5.9) 


h2 
om eg] 1— Satyr 0079], (5.10) 


and V = ve) _ A 1)r3+- ovr)| ; (5.11) 
We see from (5.10) and (5.11) that both c and V have a maximum value 
of c, when r, = 0, corresponding to p = 0. When hr, approaches an upper 
limit (m—4)z, the \alue of r, becomes very large, corresponding to short 
waves of high frequency. For large values of r, we obtain from (5.8) the 
moet Sit © = ef 1+ (y—1)r3/2yr? + O74), (5.12) 
while from the expressions (5.9) and (5.12) 
V = ef 1—(y—1)r9/2yr2+ O77 9)]. (5.13) 
Remembering our assumption that y > 1, we see that as p tends to 
infinity, c tends to the limit c, from above, but V tends to the limit c, 
from below. In any one mode, ¢ and V are continuous functions of p, 
so that V must possess a minimum which is less than c,. Denoting by 
V,, the minimum group velocity in the mth mode, we see that the 
group velocity must lie between c, and V,,. Thus for any mode the general 
shape of the dispersion curves will be as in Fig. 3, which illustrates the 
behaviour of c and V in the first mode with the constants y and 7 having 
the value 1} and 1 respectively. The existence of a minimum group 
velocity was first demonstrated by Jeffreys (1, 3). 
By inspection of (5.9), 
V > cle > c3/c,. (5.14) 
As m becomes large it can be shown from (5.9) that V tends to the limit 


c3/c except in the immediate region of r, = 0. Thus, in the higher modes 
it is to be expected that the minimum group velocity is only a little more 
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than c3/c,, and occurs near r, = (m—1)z/h. In addition, if, in (5.9) 7 is 
made large, then V will again approach c3/c, whilst if y increases, so does 
c,, as will c¢ for some particular value of r, and r,. Consequently it is to 
be expected that V,, will decrease when m, y, or 7 increase. However, in 
every case, c3/c, will be a lower bound for V. 

In the case y = 1}, c3/c, = 0-316c,. Taking + = 1, it was found that 
V,, was 0-986c,, 0-94c,, 0-91cy, and 0-84c, for m = 1, 2, 3, and 10 respec- 


tively. If y = 5, V, was found to be 0-96c,, 0-98c,, 0-86c,, and 0-85c, 
for r = 1, 2, 3, and 5 respectively. The computations were done with the 
aid of an ‘autocode’ programme for the Manchester University ‘Pegasus’ 


computer and bear out the deduction that V,, decreases for increases in 


m, y, Or T. 
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LOVE WAVES IN HYPOELASTIC BODY OF 
GRADE ZERO 
By GUNADHAR PARIA (Kharagpur, India) 
[Received 31 October 1957] 


SUMMARY 
The propagation of Love waves in the hypoelastic body of grade zero has been 
considered. It is found that the wave velocity retains the same character as that 
in the classical linearly elastic body. The difference lies in the fact that one of the 
normal stresses is not identically zero in the hypoelastic body considered, whereas 
all the normal stresses are zero in the classical elastic body. 


1. Introduction 

Tue theory of elastic response based on time rates, called hypoelasticity 
by Truesdell (1, 2, 3), has been illustrated by him and Green (4) in some 
simple problems. Thus the problems of simple extension, pure shear, 
hydrostatic pressure, and others reveal interesting phenomena charac- 
teristic of the theory. Green (5, 6) has also studied the theory for both 
incompressible and compressible bodies with reference to the theory of 
plasticity. In this paper, the propagation of Love waves in a particular 
type of hypoelasticity, that of grade zero, has been examined. It is found 
that the wave velocity possesses a feature similar to that in the classical 
linearly elastic body. There is, however, a difference; whereas in the 
classical elasticity all the normal stress components are identically zero, 
one of the normal components in hypoelasticity considered does not vanish 
identically. 


2. Fundamental equations 


The isotropic hypoelastic body of grade zero is defined by the constitu- 
tive equations (1, 2, 3, 4) 


~s v i t 
gy = 5, Ma (1) 


where ge — oe 
ot 

di; = $0; 5+;,). (3) 

In the above equations g, is the covariant metric tensor of a general 
fixed system of coordinates, v, is the covariant velocity vector, 2ysj is 
the stress tensor (« being the rigidity modulus of classical elasticity), v is 
the classical Poisson ratio, t denotes time, and a comma denotes covariant 


+ umgik — simyk, — gm*y* + attdm, (2) 


(Quart. Journ. Mech. and Applied Math., Vol. XI, Pt. 4, 1958] 





510 G. PARIA 
differentiation. The other differential equations of the theory are the 
equations of motion and continuity written in the forms 


ij 


» 
=ps 


a 
P+ (pv'), = 0, 
ot 
where p is the density and body forces are taken as zero. 


3. Love waves 
We consider a semi-infinite medium bounded by the plane z = 0, and 
the positive direction of the z-axis is taken into the medium. A layer of 
thickness 7’ of different material is spread over the surface z = 0 so that 
the upper surface of the layer is given by z = —T. Let (u,p) be the 
rigidity modulus and density respectively of the layer and (xu’,p’) be the 
corresponding values for the material below. We consider the possibility 
of the propagation of a purely transverse wave (of Love type) in the 
medium such that the disturbance penetrates only a little distance into 
the interior. Let the wave be propagated parallel to the x-axis with 
velocity c. We assume that the velocity components are 
py! vs 0, (6) 
ve vy, V(zjexp{ik(a—ct)}, (7) 
where V is a function of z only and & measures the wavelength, and where 
i= <(—1). On the right-hand side of the equation (7), only the real part 
is to be taken. The equation (5) is satisfied if 
p — constant. (8) 
Then using the relations (2) and (3) equation (1) reduces to 
0. O8 Lp 


ot 





ct cy 


where we have written s,, Srrs Sie Sryreee, CLC. 


The first three equations of (9) are satisfied if 


d & ‘ 
Sir ze Sirs ( », 
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The solution of the last two equations of (9), when relations (10) are used, 
are 


ty ‘ 


.. - din V exp{ik(a—et)} 
C (11) 


: = exp{ik(a—ct)} 


“we Dhe d 
The fourth equation of (9) is then satisfied if 
1f.. 1 /dV\* i , 
San sal! 2_ ala) Jexp(2ik(a—en)} (12) 
Two of the equations (4) are identically satisfied. The other equation 
dV c 
aes <a?) pee 13 
dr® (i (13) 
where Bp? e (14) 
P 
The solution of (13) can be written as 
V = Acos(sz)+Bsin(sz) (—T <z < 0), (15) 
V = Aexp(—s’z) (0<2z< o), (16) 


gives 


where s’ is a positive real constant. The condition that the displacement 
is continuous across z = 0 has been utilized in writing the equations (15) 
and (16). Equations (13) and (15) are satisfied if 


k 
p r 


9 


s* 


(c?—p?). (17) 
In (13) replacing B® by 8’ where 


pe 
P 
and then using (16) we have 
a= 4 K pi? ~c?)t, (19) 
B 
Since s’ is assumed to be real and positive, 
Bp’ > c. (20) 
The condition that the stress 2us,, should be continuous across z = 0 
gives, by (11), (15), and (16), 
—sBu = 8’ Ap’. 
Again the condition that 2us,, == 0 atz = —T' gives 


A sin(sT')+ Beos(sT’) = 0. 
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From (21) and (22), sutan(sT’) = s’p’. (23) 
Since the right-hand side of (23) is real and positive, it is concluded from 
this and (17) that c>B. (24) 
Thus from (20) and (24), B<ec<e. (25) 
The result (25) is well known in the classical theory of linear elasticity 
(7, 8). From (15), (21), and (23), 

V = Acos{s(z+T)}/cos(sT) (—T <2z< 0), (26) 
where v, = A at x = z= 0,t= 0. From (12), (16), and (26), we have, 
using (17) and (19), 

A*ufc?B? cos*(sT’)}-1{ B?—c? sin*{s(z+ T)}exp{2ik(a—ct)} | 

(—T <z< 0), (27) 
A2y'B’-2exp{ik(a—ct)—2p’z} (0 < z < a). (28) 
From (27) and (28) it is seen that one of the normal stresses is not 
identically zero for the hypoelastic body of grade zero, whereas in the 


o2 
8 
os vy 


9,,’e 
<F Syy 


classical linear elasticity all the normal stresses are zero identically. 
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